A tridiagonal transformation method of symmetric strip matrix based on GPU parallel acceleration

Abstract

The present invention provides a tri-diagonal transformation method of a symmetric strip matrix based on GPU parallel acceleration, which belongs to the field of high-performance computing technology. The method creates multiple thread blocks based on GPU, and the threads in the thread blocks process bump chasing in parallel; at the same time, global data storage is created on the GPU to store the striped data of the matrix to be transformed, which not only reduces the data storage volume and increases the data locality, but also is used for data synchronization between thread blocks; after all thread blocks have completed their own bump chasing, the transformed data in the global array are copied back to the corresponding positions of the main diagonal, the previous diagonal and the next diagonal of the matrix A, respectively, to obtain three striped matrices, and complete the entire bump chasing process. This method realizes bump chasing on the GPU for the first time, and realizes parallel acceleration. Compared with the existing tri-diagonalization method of symmetric strip matrices, it can greatly improve the GPU processing performance.

Description

Technical Field

[0001] The present invention belongs to the technical field of high performance computing, and in particular relates to a tridiagonal transformation method of a symmetric strip matrix based on GPU parallel acceleration. Background Art

[0002] Symmetric matrix eigenvalue decomposition is an important topic in high-performance data computing, involving multiple fields such as artificial intelligence and scientific computing. The commonly used method for symmetric matrix eigenvalue decomposition is a three-stage decomposition algorithm, which is divided into the symmetric matrix striping stage, the symmetric striped matrix tridiagonalization stage and the eigenvalue calculation stage, which respectively completes the three tasks of striping decomposition of symmetric matrices, decomposing striped symmetric matrices into tridiagonalization, and using the QR iterative algorithm or divide-and-conquer method to calculate the final eigenvalue. The performance bottleneck of the existing algorithm is reflected in the first two stages, which account for more than 95% of the time in the commonly used algorithm package CUDA cuSolver.

[0003] The commonly used method for tridiagonal decomposition of symmetric striped matrices is bulge chasing. However, this method is generally considered to be a memory bottleneck. It is believed that the bottleneck of the algorithm mainly exists in memory transfer, which is a memory bottleneck problem rather than a computing bottleneck problem. Therefore, there is no version implemented on the GPU so far. Summary of the invention

[0004] With the continuous development of hardware acceleration technology, NVIDIA GPU has achieved great success in the fields of artificial intelligence, parallel computing, etc. Its SIMD (Single Instruction / Multiple Data) architecture can run multiple thread blocks at the same time, and each thread block can run multiple threads at the same time. Such extremely high parallelism provides huge potential for the parallel optimization of multiple algorithms. Although the bump chasing algorithm can only perform serial calculations in a single bump chasing process, its multi-bump chasing process has huge parallel potential. The present invention realizes a novel and efficient GPU parallel algorithm, and implements the entire bump chasing process in parallel. Comparing MAGMA (an open source linear algebra software library from the University of Tennessee, USA) and cuSolver (a linear algebra software library released by NVIDIA), the performance is improved by 7.5 times and 10.1 times respectively.

[0005] In order to solve the above technical problems, the specific technical solution of the tridiagonal transformation method of a symmetric strip matrix based on GPU parallel acceleration of the present invention is as follows:

[0006] A symmetric strip matrix tridiagonal transformation method based on GPU parallel acceleration, for an n*n symmetric strip matrix A, whose strip width is b, performs a tridiagonal transformation, comprising the following steps:

[0007] Step 1: Create a global array dS[2*b*n] of size 2*b*n on the GPU for data synchronization among thread blocks in step 2; set variables i and j to represent the number of rows and columns of elements in the matrix A to be copied, and copy the elements in the matrix A that satisfy j≤i≤j+b to the array dS;

[0008] Specifically, initialize j = 0. The specific steps of copying the elements in matrix A that satisfy j≤i≤j+b to array dS are as follows:

[0009] Step 1.1: Initialize i=0, set the number of data to be copied count=min(n,j+b+1)-j, and jump to step 1.2;

[0010] Step 1.2: Determine whether i is less than count. If so, set dS[i+j*2*b]=A[i+j+j*n], i=i+1, and jump to step 1.2. Otherwise, jump to step 1.3.

[0011] Step 1.3: Let j = j + 1, determine whether j is less than n, if so, jump to step 1.1, otherwise jump to step 2.

[0012] Step 2: Create n-2 thread blocks on the GPU. Each thread block processes the bump chase from the 1st column to the n-2th column. Each thread block includes k*k threads. The k*k threads process a single bump chase together. Initialize the array com[n-2] for synchronizing the processing positions of each thread block and for recording the current position opRow of each thread block for bump chasing. When initialized, com[n-2]={0}.

[0013] Step 3: Each thread block processes the corresponding bump chase separately; let the current position processed by each thread block be opRow, and the number of each thread block be bIndx. The Householder transformation process of each bump chase involves three parts, namely B1, S and B2; the number of rows and columns of B1 are b1Row and b1Col respectively, and the number of rows and columns of S are sRow and sCol respectively. The number of rows and columns of B2 are b2Row and b2Col respectively; specifically, in a bump chase, B1 is the part of the matrix A with rows from opRow to opRow+b1Row-1 and columns from opRow-b1Col to opRow-1; S is the part of the matrix A with rows from opRow to opRow+sRow-1 and columns from opRow to opRow+sCol-1; B2 is the part of the matrix A with rows from opRow+sRow to opRow+sRow+b2Row-1 and columns from opRow to opRow+b2Col-1; in each thread block, create shared memory arrays S_B1[b*b], S_S[b*b], and S_B2[b*b] to store B1, S, and B2 respectively; let b2Col=1, opRow=bIndx+1; the specific process of bump chase is as follows:

[0014] Step 3.1: If the current thread block processing event meets the following conditions:

[0015] bIndx≠0

[0016] And opRow+2*b>com[bInx-1]

[0017] Then loop judgment is performed, otherwise jump to 3.2; where com[bInx-1] represents the bInxth element in the array;

[0018] Step 3.2:

[0019] b1Col=b2Col

[0020] b1Row=min(b,n-opRow)

[0021] b2Col=sCol =sRow= b1Row

[0022] b2Row=min(b, n-sRow-opRow)

[0023] Copy the corresponding parts of array dS to S_B1, S_S and S_B2, and jump to step 3.3;

[0024] Step 3.3: First, perform Householder transformation on the first column of S_B1 to obtain the corresponding Householder vector u.

[0025] Transform S_B1, S_S, and S_B2 respectively, where I is the unit matrix with a size of b1Row *b1Row; jump to step 3.4;

[0026] Step 3.4: Copy S_B1, S_S and S_B2 to the corresponding parts of array dS, and jump to step 3.5;

[0027] Step 3.5: Order

[0028] opRow=opRow+b

[0029] com[bInx] = opRow

[0030] Determine whether b2Row is less than 2. If so, go to step 4. If not, jump to step 3.1; where com[bInx] represents the bInx+1th element in the array.

[0031] Step 4: Wait until all thread blocks have completed their own bump chasing, copy the first row in array dS to the main diagonal of matrix A, copy the second row in array dS to the upper and lower diagonals of matrix A to obtain three striped matrices, and complete the entire bump chasing process.

[0032] Specifically, the step 3.2 is as follows:

[0033] Step 3.2.1: Set variables dS_B1, dS_S and dS_B2 to the corresponding positions of B1, S and B2 in dS to be copied respectively; let dS_B1=dS+b1Col + (opRow - b1Col)*2*b, dS_S=dS+ opRow*2*b, dS_B2=dS+sRow+opRow*2*b; let i and j be the number of rows and columns of the data to be copied respectively, let j=0, and jump to step 3.2.2;

[0034] Step 3.2.2: Determine whether j is less than b1Col. If so, set i=0 and jump to 3.2.3; otherwise, set j=0 and jump to 3.2.4;

[0035] Step 3.2.3: Determine whether i is less than b1Row. If so, set S_B1[i+j*b]= dS_B1[i+(ij)*b], i++, and jump to 3.2.3; otherwise, set j=j+1 and jump to 3.2.2;

[0036] Step 3.2.4: Determine whether j is less than sCol. If so, set i=0 and jump to 3.2.5; otherwise, set j=0 and jump to 3.2.6;

[0037] Step 3.2.5: Determine whether i is less than sRow. If so, set S_S[i+j*b]= dS_S[i+(ij)*b], i++, and jump to 3.2.5; otherwise, set j=j+1 and jump to 3.2.4;

[0038] Step 3.2.6: Determine whether j is less than b2Col. If so, set i=0 and jump to 3.2.7; otherwise, complete step 3.2 and jump to step 3.3;

[0039] Step 3.2.7: Determine whether i is less than b2Row. If so, set S_B2[i+j*b]= dS_B2[i+(ij)*b], i++, and jump to 3.2.7; otherwise, set j=j+1 and jump to 3.2.6;

[0040] Specifically, step 3.4 is exactly opposite to step 3.2, as follows:

[0041] Step 3.4.1: Set variables dS_B1, dS_S and dS_B2 to the corresponding positions of B1, S and B2 in dS to be copied; let dS_B1=dS+b1Col + (opRow - b1Col) *2*b, dS_S=dS+ opRow*2*b, dS_B2=dS+sRow+opRow*2*b; let i and j be the number of rows and columns of the data to be copied, let j=0, and jump to step 3.4.2;

[0042] Step 3.4.2: Determine whether j is less than b1Col. If so, set i=0 and jump to 3.4.3; otherwise, set j=0 and jump to 3.4.4;

[0043] Step 3.4.3: Determine whether i is less than b1Row. If so, set dS_B1[i+(ij)*b]=S_B1[i+j*b], i++, and jump to 3.4.3; otherwise, set j=j+1 and jump to 3.4.2;

[0044] Step 3.4.4: Determine whether j is less than sCol. If so, set i=0 and jump to 3.4.5; otherwise, set j=0 and jump to 3.4.6;

[0045] Step 3.4.5: Determine whether i is less than sRow. If so, set dS_S[i+(ij)*b]=S_S[i+j*b], i++, and jump to 3.4.5; otherwise, set j=j+1 and jump to 3.4.4;

[0046] Step 3.4.6: Determine whether j is less than b2Col. If so, set i = 0 and jump to 3.4.7; otherwise, complete Step 3.4 and jump to Step 3.5.

[0047] Step 3.4.7: Determine whether i is less than b2Row. If so, set dS_B2[i+(i-j)*b]=S_B2[i+j*b], increment i, and jump to 3.4.7; otherwise, set j = j + 1 and jump to 3.4.6.

[0048] Specifically, Step 3.3 is as follows:

[0049] Step 3.3.1: Divide the k * k threads created in each thread block into k warps, numbered p (0 ≤ p < k); each warp contains k threads, numbered q (0 ≤ q < k); jump to 3.3.2.

[0050] Step 3.3.2: Let the variable col be the column number of S_B1 to be subjected to the Householder transformation, and set col = p; jump to Step 3.3.3.

[0051] Step 3.3.3: Determine whether col is less than b1Col. If so, the k threads in warp p perform the Householder transformation on the col column of S_B1, set col = col + k, and jump to 3.3.3; otherwise, jump to Step 3.3.4.

[0052] Step 3.3.4: Let the variable col be the column number of S_S to be subjected to the Householder transformation, and set col = p; jump to Step 3.3.5.

[0053] Step 3.3.5: Determine whether col is less than sCol. If so, the k threads in warp p perform the Householder transformation on the col column of S_S, set col = col + k, and jump to 3.3.5; otherwise, jump to Step 3.3.6.

[0054] Step 3.3.6: Let the variable row be the row number of S_S to be subjected to the Householder transformation, and set row = p; jump to Step 3.3.7.

[0055] Step 3.3.7: Determine whether row is less than b1Row. If so, the k threads in warp p perform the Householder transformation on the row row of S_S, set row = row + k, and jump to 3.3.7; otherwise, jump to Step 3.3.8.

[0056] Step 3.3.8: Let the variable col be the number of columns of S_B2 to be Householder transformed, and let col=p; jump to step 3.3.9;

[0057] Step 3.3.9: Determine whether col is less than b2Col. If so, the k threads in the thread bundle p perform Householder transformation on the col column of S_B2, set col=col+k, and jump to 3.3.9; otherwise, jump to step 3.4;

[0058] Specifically, the step 4 is as follows:

[0059] Step 4.1: Let i = 0, which is the number of columns of the array dS to be copied;

[0060] Step 4.2: Determine whether i is less than n. If so, jump to step 4.3; otherwise, end;

[0061] Step 4.3: Determine whether i is equal to 0 and equal to n-1. If not, let A[i+i*n]=dS[i*2*b], A[i+1+i*n]=dS[1+i*2*b], A[i-1+i*n]=dS[1+i*2*b], i=i+1 jump to step 4.2; Determine whether i is equal to n-1. If so, let A[i+i*n]=dS[i*2*b], A[i-1+i*n]=dS[1+i*2*b], i=i+1, jump to step 4.2; Determine whether i is equal to 0. If so, let A[i+i*n]=dS[i*2*b], A[i+1+i*n]=dS[1+i*2*b], i=i+1, jump to step 4.2.

[0062] The beneficial effects of the present invention are as follows:

[0063] The present invention provides a symmetric strip matrix tridiagonalization transformation method based on GPU parallel acceleration. This method realizes bump chasing on GPU for the first time and realizes parallel acceleration. Compared with the existing symmetric strip matrix tridiagonalization method, it can greatly improve GPU processing performance and has high computational efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 It is a schematic diagram of the process of the present invention;

[0065] Figure 2 It is a schematic diagram of the original symmetric striped matrix;

[0066] Figure 3 Schematic diagram of the orthogonal transformation matrix for the first bump chase on the first column of the matrix;

[0067] Figure 4This is a schematic diagram of performing Householder transformation on the 2nd to n-2nd columns of a matrix;

[0068] Figure 5 Schematic diagram of the results of bump chasing for columns 2 to n-2 of the matrix;

[0069] Figure 6 It is a three-strip transformation diagram of the symmetric matrix A;

[0070] Figure 7 Schematic diagram of the Householder transformation process in the method of the present invention. DETAILED DESCRIPTION

[0071] In order to better understand the purpose, structure and function of the present invention, the following is a further detailed description of an efficient bump chasing method using GPU parallel acceleration of the present invention in conjunction with the accompanying drawings.

[0072] First, we take a symmetric striped matrix with a matrix size of n=15 and a striping width of b=3 as an example to illustrate the convex block chasing method. The original symmetric striped matrix A is as follows: Figure 2 As shown, each square in the figure codes 1 element of the matrix A, the blank square represents the 0 element, and x represents any non-zero element.

[0073] A(i:j,h:k) is used to extract the elements from the i-th row to the j-th row and the h-th column to the k-th column in the matrix A. First, the first bump chasing process is performed on the first column of the matrix, and the non-zero elements below the second diagonal are eliminated. The three elements of A(1:1,2:4) are used as the vectors to be transformed, and the Householder transformation is performed to eliminate the two elements of A(1:1,3:4). The Householder vector u obtained by the Householder transformation is used to perform a similarity transformation on A:

[0074] It will generate non-zero elements at A(2:2,6:7) and A(3:3,7:7), and these three elements will form a convex block (it will also generate corresponding convex blocks in the upper triangular part of the matrix A, but no repeated calculations are performed due to symmetry), such as Figure 3 shown.

[0075] Then, the three elements of A(2:2,5:7) are used as the vectors to be transformed, and the Householder transformation is performed on them, which eliminates the two elements of A(2:2,6:7). At the same time, it generates three new non-zero elements at A(5:5,8:9) and A(6:6,9:9), which is equivalent to generating new convex blocks, such as Figure 4 The above process is repeated until all bumps are chased out of the entire matrix.

[0076] Similarly, perform convex block chasing on the 2nd to n-2nd columns of matrix A. Eliminate the non-zero elements below the diagonal of the 2nd to n-2nd columns respectively, and get Figure 5 In the matrix shown in FIG. 1 , the blank part is a 0 element, the white part is a non-zero element, and the gray part is an element not involved in the calculation; using symmetry, the lower triangular part of the symmetric matrix A is copied to the upper triangular part of the symmetric matrix A to obtain a three-strip transformation form of the symmetric matrix A, as shown in FIG. Figure 6 As shown, the entire bump chasing process is completed.

[0077] The following is a specific embodiment of the present invention, a symmetric strip matrix tridiagonal transformation method based on GPU parallel acceleration, for a 32768*32768 symmetric strip matrix A, whose strip width b is 32, a tridiagonal transformation is performed, such as Figure 1 As shown, the following steps are included:

[0078] Step 1: Create a global array dS[2*32*32768] of size 2*32*32768 on the GPU for data synchronization of each thread block in step 2; set variables i and j to represent the number of rows and columns of the elements in the matrix A to be copied, and copy the elements in the matrix A that satisfy j≤i≤j+32 to the array dS;

[0079] Specifically, initialize j = 0. The specific steps of copying the elements in matrix A that satisfy j≤i≤j+32 to array dS are as follows:

[0080] Step 1.1: Initialize i=0, set the number of data to be copied count=min(n,j+32+1)-j, and jump to step 1.2;

[0081] Step 1.2: Determine whether i is less than count. If so, set dS[i+j*2*32]=A[i+j+j*32768], i=i+1, and jump to step 1.2. Otherwise, jump to step 1.3.

[0082] Step 1.3: Let j = j + 1, and determine whether j is less than 32768. If so, jump to step 1.1, otherwise jump to step 2.

[0083] Step 2: Perform basic initialization. Create 32768-2=32766 thread blocks on the GPU. Each thread block processes bump chasing from the 1st column to the 32768-2nd column. Each thread block includes 32*32 threads. 32*32 threads process a single bump chasing together. Initialize the array com[32768-2] for each thread block to synchronize the processing position and to record the current position opRow of each thread block for bump chasing. During initialization, com[32768-2]={0}.

[0084] Step 2: Each thread block processes its corresponding convex block chase separately; let the current position processed by each thread block be opRow, and the Householder transformation process of each convex block chase involves three parts, namely B1, S and B2, as shown in Figure 7 As shown, B1 is the orange part, S is the green part, B2 is the pink part, the blank part is the 0 element, the white part is the non-zero element, and the gray part is the element not involved in the calculation. Let the number of rows and columns of B1 be b1Row and b1Col respectively, let the number of rows and columns of S be sRow and sCol respectively, let the number of rows and columns of B2 be b2Row and b2Col respectively; specifically, in a bump chase, B1 is the part of the matrix A with rows from opRow to opRow+b1Row-1 and columns from opRow-b1Col to opRow-1; S is the part of the matrix A with rows from opRow to opRow+sRow-1 and columns from opRow to opRow+sCol-1; B2 is the part of the matrix A with rows from opRow+sRow to opRow+sRow+b2Row-1 and columns from opRow to opRow+b2Col-1.

[0085] In each thread block, create shared memory arrays S_B1[b*b], S_S[b*b], and S_B2[b*b] to store B1, S, and B2 respectively; let b2Col=1, opRow=bIndx+1; the specific process of convex chasing is as follows:

[0086] Step 3.1: If the current thread block processing event meets the following conditions:

[0087] bIndx≠0

[0088] And opRow+2*32>com[bInx-1]

[0089] Then loop judgment is performed, otherwise jump to 3.2; where com[bInx-1] represents the bInxth element in the array, that is, the processing position of the bInxth bump chase. The bInx+1th bump chase will poll the processing position of the bInxth bump chase (that is, com[bInx-1]). When the data dependency is met, the bInx+1th bump chase can be processed.

[0090] Step 3.2: Order

[0091] b1Col=b2Col

[0092] b1Row=min(32,32768-opRow)

[0093] b2Col=sCol =sRow= b1Row

[0094] b2Row=min(32,32768-sRow-opRow)

[0095] Copy the corresponding part of array dS to S_B1, S_S and S_B2, and jump to step 3.3; the specific memory update is as follows:

[0096] Step 3.2.1: Set variables dS_B1, dS_S and dS_B2 to the corresponding positions of B1, S and B2 in dS to be copied respectively; let dS_B1=dS+b1Col + (opRow - b1Col)*2*32, dS_S=dS+ opRow*2*32, dS_B2=dS+sRow+opRow*2*32; let i and j be the number of rows and columns of the data to be copied respectively, let j=0, jump to step 3.2.2;

[0097] Step 3.2.2: Determine whether j is less than b1Col. If so, set i=0 and jump to 3.2.3; otherwise, set j=0 and jump to 3.2.4;

[0098] Step 3.2.3: Determine whether i is less than b1Row. If so, set S_B1[i+j*32]= dS_B1[i+(ij)*32], i++, and jump to 3.2.3; otherwise, set j=j+1 and jump to 3.2.2;

[0099] Step 3.2.4: Determine whether j is less than sCol. If so, set i=0 and jump to 3.2.5; otherwise, set j=0 and jump to 3.2.6;

[0100] Step 3.2.5: Determine whether i is less than sRow. If so, let S_S[i + j * 32] = dS_S[i + (i - j) * 32], increment i by 1, and jump to 3.2.5; otherwise, let j = j + 1 and jump to 3.2.4;

[0101] Step 3.2.6: Determine whether j is less than b2Col. If so, let i = 0 and jump to 3.2.7; otherwise, complete Step 3.2 and jump to Step 3.3;

[0102] Step 3.2.7: Determine whether i is less than b2Row. If so, let S_B2[i + j * 32] = dS_B2[i + (i - j) * 32], increment i by 1, and jump to 3.2.7; otherwise, let j = j + 1 and jump to 3.2.6;

[0103] Step 3.3: First, perform a Householder transformation on the first column of S_B1 to obtain the corresponding Householder vector u, and let

[0104] Perform transformations on the three parts of S_B1, S_S, and S_B2 respectively. Here, I is the identity matrix with size b1Row * b1Row; the specific steps are as follows:

[0105] Step 3.3.1: Divide the k * k threads created in each thread block into k thread warps, numbered p (0 ≤ p < k); each thread warp contains k threads, numbered q (0 ≤ q < k); jump to 3.3.2;

[0106] Step 3.3.2: Let the variable col be the column number of S_B1 to be subjected to the Householder transformation, and let col = p; jump to Step 3.3.3;

[0107] Step 3.3.3: Determine whether col is less than b1Col. If so, the k threads in thread warp p perform a Householder transformation on the col column of S_B1, let col = col + k, and jump to 3.3.3; otherwise, jump to Step 3.3.4;

[0108] Step 3.3.4: Let the variable col be the column number of S_S to be subjected to the Householder transformation, and let col = p; jump to Step 3.3.5;

[0109] Step 3.3.5: Determine whether col is less than sCol. If so, the k threads in thread warp p perform a Householder transformation on the col column of S_S, let col = col + k, and jump to 3.3.5; otherwise, jump to Step 3.3.6;

[0110] Step 3.3.6: Let the variable row be the number of rows of S_S to be Householder transformed, and let row = p; jump to step 3.3.7;

[0111] Step 3.3.7: Determine whether row is less than b1Row. If so, the k threads in thread bundle p perform Householder transformation on the row of S_S, set row=row+k, and jump to 3.3.7; otherwise, jump to step 3.3.8;

[0112] Step 3.3.8: Let the variable col be the number of columns of S_B2 to be Householder transformed, and let col=p; jump to step 3.3.9;

[0113] Step 3.3.9: Determine whether col is less than b2Col. If so, the k threads in the thread bundle p perform Householder transformation on the col column of S_B2, set col=col+k, and jump to 3.3.9; otherwise, jump to step 3.4;

[0114] Step 3.4: Copy S_B1, S_S and S_B2 to the corresponding parts of the array dS; specifically, the step 3.4 is exactly the opposite of step 3.2, as follows:

[0115] Step 3.4.1: Set the variables dS_B1, dS_S and dS_B2 to the corresponding positions of B1, S and B2 in dS to be copied respectively; let dS_B1=dS+b1Col + (opRow - b1Col) *2*32, dS_S=dS+ opRow*2*32, dS_B2=dS+sRow+opRow*2*32; let i and j be the number of rows and columns of the data to be copied respectively, let j=0, jump to step 3.4.2;

[0116] Step 3.4.2: Determine whether j is less than b1Col. If so, set i=0 and jump to 3.4.3; otherwise, set j=0 and jump to 3.4.4;

[0117] Step 3.4.3: Determine whether i is less than b1Row. If so, set dS_B1[i+(ij)*32]=S_B1[i+j*32], i++, and jump to 3.4.3; otherwise, set j=j+1 and jump to 3.4.2;

[0118] Step 3.4.4: Determine whether j is less than sCol. If so, set i=0 and jump to 3.4.5; otherwise, set j=0 and jump to 3.4.6;

[0119] Step 3.4.5: Determine whether i is less than sRow. If so, set dS_S[i+(ij)*32]=S_S[i+j*32], i++, and jump to 3.4.5; otherwise, set j=j+1 and jump to 3.4.4;

[0120] Step 3.4.6: Determine whether j is less than b2Col. If so, set i=0 and jump to 3.4.7. Otherwise, complete step 3.4 and jump to step 3.5.

[0121] Step 3.4.7: Determine whether i is less than b1Row. If so, set dS_B1[i+(ij)*32]=S_B1[i+j*32], i++, and jump to 3.4.7; otherwise, set j=j+1 and jump to 3.4.6;

[0122] Step 3.5: Order

[0123] opRow=opRow+32

[0124] com[bInx] = opRow

[0125] Determine whether b2Row is less than 2. If so, go to step 4. If not, jump to step 3.1; where com[bInx] represents the bInx+1th element in the array.

[0126] Step 4: Wait until all thread blocks have completed their own bump chasing, copy the first row in array dS to the main diagonal of matrix A, copy the second row in array dS to the upper and lower diagonals of matrix A to obtain three striped matrices, and complete the entire bump chasing process; the specific steps are as follows:

[0127] Step 4.1: Let i = 0, which is the number of columns of the array dS to be copied;

[0128] Step 4.2: Determine whether i is less than 32768. If so, jump to step 4.3; otherwise, end;

[0129] Step 4.3: Determine whether i is equal to 0 and equal to 32768-1. If not, set A[i+i*32768]=dS[i*2*32], A[i+1+i*32768]=dS[1+i*2*32], A[i-1+i*32768]=dS[1+i*2*32], i=i+1 and jump to step 4.2; Determine whether i is equal to 32768-1. If so, set A[ i+i*32768]=dS[i*2*32], A[i-1+i*32768]=dS[1+i*2*32], i=i+1, jump to step 4.2; determine whether i is equal to 0. If A[i+i*32768]=dS[i*2*32], A[i+1+i*32768]=dS[1+i*2*32], i=i+1, jump to step 4.2.

[0130] The commonly used method for symmetric matrix eigenvalue decomposition is a three-stage decomposition algorithm, which is divided into a symmetric matrix striping stage, a symmetric striping matrix tridiagonalization stage, and an eigenvalue calculation stage. The method of the present invention is used to perform tridiagonalization of symmetric striping matrices, which will greatly improve the performance of bump chasing. On the NVIDIA H100 GPU, each block contains 32*32 threads. For the symmetric striping matrix b=32, the performance of the MAGMA single bump chasing stage and the first two stages of CUDAcuSolver (cuSolver does not provide an interface for bump chasing alone, so only the first two stages can be used for comparison) are improved by 7.5 times and 10.1 times respectively. As shown in Table 1:

[0131] Table 1

[0132]

[0133] It is to be understood that the present invention is described by some embodiments, and it is known to those skilled in the art that various changes or equivalent substitutions may be made to these features and embodiments without departing from the spirit and scope of the present invention. In addition, under the teachings of the present invention, these features and embodiments may be modified to adapt to specific circumstances and materials without departing from the spirit and scope of the present invention. Therefore, the present invention is not limited by the specific embodiments disclosed herein, and all embodiments falling within the scope of the claims of this application are within the scope of protection of the present invention.

Claims

1. A symmetric strip matrix tridiagonal transformation method based on GPU parallel acceleration, for an n*n symmetric strip matrix A, whose strip width is b, performs a tridiagonal transformation, characterized in that: The following steps are involved: Step 1: Create a global array dS[2*b*n] of size 2*b*n on the GPU for data synchronization among thread blocks in step 2; Set variables i and j to represent the number of rows and columns of the elements in the matrix A to be copied, and copy the elements in the matrix A that satisfy j≤i≤j+b to the array dS; Step 2: Create n-2 thread blocks on the GPU. Each thread block processes the bump chase from the 1st column to the n-2th column. Each thread block includes k*k threads. The k*k threads process a single bump chase together. Initialize the array com[n-2] for synchronizing the processing positions of each thread block and for recording the current position opRow of each thread block for bump chasing. When initialized, com[n-2] = {0}. Step 3: Each thread block processes the corresponding convex block chase separately; Let the current position processed by each thread block be opRow, and the number of each thread block be bIndx. The Householder transformation process of each bump chase involves three parts, namely B1, S and B2; the number of rows and columns of B1 are b1Row and b1Col respectively, the number of rows and columns of S are sRow and sCol respectively, and the number of rows and columns of B2 are b2Row and b2Col respectively; specifically, in a bump chase, B1 is the part of the matrix A with rows from opRow to opRow+b1Row-1 and columns from opRow-b1Col to opRow-1; S is the portion of the matrix A with rows from opRow to opRow+sRow-1 and columns from opRow to opRow+sCol-1; B2 is the part of the matrix A with rows from opRow+sRow to opRow+sRow+b2Row-1 and columns from opRow to opRow+b2Col-1. In each thread block, create shared memory arrays S_B1[b*b], S_S[b*b], and S_B2[b*b] to store B1, S, and B2 respectively. Set b2Col=1, opRow=bIndx+1. The specific process of convex chasing is as follows: Step 3.1: If the current thread block processing event meets the following conditions: bIndx≠0 And opRow+2*b>com[bInx-1] Then loop judgment is performed, otherwise jump to 3.2; where com[bInx-1] represents the bInxth element in the array; Step 3.2: Order b1Col=b2Col b1Row=min(b,n-opRow) b2Col=sCol=sRow=b1Row b2Row=min(b,n-sRow-opRow) Copy the corresponding part of the array dS to S_B1, S_S and S_B2, and jump to step 3.3; specifically, the step 3.2 is as follows: Step 3.2.1: Set variables dS_B1, dS_S, and dS_B2 to the corresponding positions of B1, S, and B2 in dS to be copied respectively; let dS_B1 = dS + b1Col+(opRow - b1Col)*2*b, dS_S = dS + opRow*2*b, dS_B2 = dS + sRow + opRow*2*b; let i and j be the number of rows and columns of the data to be copied respectively, let j = 0, and jump to Step 3.2.2; Step 3.2.2: Determine whether j is less than b1Col. If so, let i = 0 and jump to 3.2.3; otherwise, let j = 0 and jump to 3.2.4; Step 3.2.3: Determine whether i is less than b1Row. If so, let S_B1[i + j*b] = dS_B1[i+(i - j)*b], i++, and jump to 3.2.3; Otherwise, let j = j + 1 and jump to 3.2.2; Step 3.2.4: Determine whether j is less than sCol. If so, let i = 0 and jump to 3.2.5; otherwise, let j = 0 and jump to 3.2.6; Step 3.2.5: Determine whether i is less than sRow. If so, let S_S[i + j*b] = dS_S[i+(i - j)*b], i++, and jump to 3.2.5; otherwise, let j = j + 1 and jump to 3.2.4; Step 3.2.6: Determine whether j is less than b2Col. If so, let i = 0 and jump to 3.2.7; Otherwise, complete Step 3.2 and jump to Step 3.3; Step 3.2.7: Determine whether i is less than b2Row. If so, let S_B2[i + j*b] = dS_B2[i+(i - j)*b], i++, and jump to 3.2.7; otherwise, let j = j + 1 and jump to 3.2.6; Step 3.3: First, perform a Householder transformation on the first column of S_B1 to obtain the corresponding Householder vector u, and let Perform transformations on the three parts of S_B1, S_S, and S_B2 respectively, where I is the identity matrix with size b1Row*b1Row; jump to Step 3.4; The specific steps of Step 3.3 are as follows: Step 3.3.1: Divide the k*k threads created in each thread block into k warps, numbered p(0 ≤ p < k); each warp contains k threads, numbered q(0 ≤ q < k); jump to 3.3.2; Step 3.3.2: Let the variable col be the column number of S_B1 to be subjected to the Householder transformation, and let col = p; jump to Step 3.3.3; Step 3.3.3: Determine whether col is less than b1Col. If so, the k threads in warp p perform a Householder transformation on the col column of S_B1, let col = col + k, and jump to 3.3.3; otherwise, jump to Step 3.3.4; Step 3.3.4: Let the variable col be the column number of S_S to be subjected to the Householder transformation, and let col = p; jump to Step 3.3.5; Step 3.3.5: Determine whether col is less than sCol. If so, the k threads in the thread bundle p perform Householder transformation on the col column of S_S, set col = col + k, and jump to 3.3.5; otherwise, jump to step 3.3.6; Step 3.3.6: Let the variable row be the number of rows of S_S to be Householder transformed, and let row = p; jump to step 3.3.7; Step 3.3.7: Determine whether row is less than b1Row. If so, the k threads in thread bundle p perform Householder transformation on row of S_S, set row = row + k, and jump to 3.3.7; otherwise, jump to step 3.3.8; Step 3.3.8: Let the variable col be the number of columns of S_B2 to be Householder transformed, and let col = p; jump to step 3.3.9; Step 3.3.9: Determine whether col is less than b2Col. If so, the k threads in the thread bundle p perform Householder transformation on the col column of S_B2, set col = col + k, and jump to 3.3.9; otherwise, jump to step 3.4; Step 3.4: Copy S_B1, S_S and S_B2 to the corresponding parts of array dS, and jump to step 3.5; step 3.4 is exactly the opposite of step 3.2, as follows: Step 3.4.1: Set the variables dS_B1, dS_S and dS_B2 to the corresponding positions of B1, S and B2 in dS to be copied; let dS_B1 ​​= dS + b1Col + (opRow-b1Col)*2*b, dS_S = dS + opRow*2*b, dS_B2=dS+sRow+opRow*2*b; let i and j be the number of rows and columns of the data to be copied, let j=0, and jump to step 3.4.2; Step 3.4.2: Determine whether j is less than b1Col. If so, set i=0 and jump to 3.4.3; otherwise, set j=0 and jump to 3.4.4; Step 3.4.3: Determine whether i is less than b1Row. If so, set dS_B1[i+(ij)*b]=S_B1[i+j*b], i++, and jump to 3.4.3; otherwise, set j=j+1, and jump to 3.4.2; Step 3.4.4: Determine whether j is less than sCol. If so, set i = 0 and jump to 3.4.5; otherwise, set j = 0 and jump to 3.4.6; Step 3.4.5: Determine whether i is less than sRow. If so, set dS_S[i+(ij)*b]=S_S[i+j*b], i++, and jump to 3.4.5; otherwise, set j=j+1, and jump to 3.4.4; Step 3.4.6: Determine whether j is less than b2Col. If so, set i=0 and jump to step 3.4.

7. Otherwise, complete step 3.4 and jump to step 3.

5. Step 3.4.7: Determine whether i is less than b2Row. If so, set dS_B2[i+(ij)*b]=S_B2[i+j*b], i++, and jump to 3.4.7; otherwise, set j=j+1, and jump to 3.4.6; Step 3.5: Order opRow=opRow+b com[bInx]=opRow Determine whether b2Row is less than 2. If so, go to step 4. If not, jump to step 3.

1. Where com[bInx] represents the bInx+1th element in the array. Step 4: Wait until all thread blocks have completed their own bump chasing, copy the first row in array dS to the main diagonal of matrix A, copy the second row in array dS to the upper and lower diagonals of matrix A to obtain three striped matrices, and complete the entire bump chasing process.

2. The method for tridiagonalization of symmetric strip matrix based on GPU parallel acceleration according to claim 1, characterized in that: The step 1 is specifically as follows: Initialize j = 0; copy the elements in matrix A that satisfy j≤i≤j+b to array dS as follows: Step 1.1: Initialize i=0, set the number of data to be copied count=min(n,j+b+1)-j, and jump to step 1.2; Step 1.2: Determine whether i is less than count. If so, set dS[i+j*2*b]=A[i+j+j*n], i=i+1, and jump to step 1.2; otherwise, jump to step 1.3; Step 1.3: Let j=j+1, determine whether j is less than n, if so, jump to step 1.1, otherwise jump to step 2.

3. The method for tridiagonalization of symmetric strip matrix based on GPU parallel acceleration according to claim 1, characterized in that: The step 4 is specifically as follows: Step 4.1: let i=0, which is the number of columns of the array dS to be copied; Step 4.2: Determine whether i is less than n. If so, jump to step 4.3; otherwise, end; Step 4.3: Determine whether i is equal to 0 and equal to n-1. If not, let A[i+i*n]=dS[i*2*b], A[i+1+i*n]=dS[1+i*2*b], A[i-1+i*n]=dS[1+i*2*b], i=i+1 and jump to step 4.2; Determine whether i is equal to n-1. If so, let A[i+i*n]=dS[i*2*b], A[i-1+i*n]=dS[1+i*2*b], i=i+1 and jump to step 4.2; Determine whether i is equal to 0. If so, let A[i+i*n]=dS[i*2*b], A[i+1+i*n]=dS[1+i*2*b], i=i+1 and jump to step 4.

2.

4. The method for tridiagonalization of symmetric strip matrix based on GPU parallel acceleration according to claim 1, characterized in that: The GPU has a single instruction multiple data SIMD architecture.

Images