A method, system and device for eigenvalue decomposition of heterogeneous symmetric matrices based on GPU parallel acceleration

Abstract

The present invention provides a method, system and device for eigenvalue decomposition of heterogeneous symmetric matrices based on GPU parallel acceleration, which belongs to the field of high-performance computing technology. The present invention breaks the design that the five sub-processes in the traditional symmetric matrix eigenvalue decomposition process can only be strictly executed in serial, solves the problem of idle CPU and GPU computing resources in serial algorithms, and greatly improves computing efficiency. At the same time, the present invention proposes a novel BC-Back method based on GPU implementation for the BC-Back stage, which transfers the BC-Back process from the CPU with limited parallel capability to the GPU that can be large-scale parallel, greatly improving the parallel efficiency of the BC-Back stage. The present invention makes full use of the characteristics of heterogeneous architecture and GPU parallelism, realizes the parallelism between the eigenvalue solution sub-processes of the symmetric matrix and the internal parallelism of BC-Back, makes full use of CPU and GPU resources, and greatly improves the computing performance of the computer for the eigenvalue solution process of the symmetric matrix.

Description

Technical Field

[0001] The present invention belongs to the field of high-performance computing technology, and in particular relates to a method, system and device for eigenvalue decomposition of heterogeneous symmetric matrices based on GPU parallel acceleration. Background Art

[0002] Symmetric matrix eigenvalue decomposition is an important topic in high-performance data computing today, involving multiple fields such as artificial intelligence and scientific computing. Symmetric matrix eigenvalue decomposition involves two versions: one that requires finding the eigenvalue vector and the other that does not. Usually in the field of artificial intelligence or high-performance computing, the version that requires finding the eigenvector has a greater demand. The present invention mainly relates to the version that requires finding the eigenvector. The currently commonly used method for symmetric matrix eigenvalue decomposition (finding the eigenvector) is a five-stage decomposition method, which is divided into 1) Successive Band Reduction (SBR) stage; 2) Bugle Chasing (BC); 3) Divide and Conquer (DC); 4) Back Transformation of Bugle Chasing (BC-Back) and 5) Back Transformation of Successive Band Reduction (SBR-Back). In view of the strict data dependency between the five stages of the existing decomposition method, all five stages are processed serially, such as Figure 1 As shown in the figure above, these five processes can only be run independently on the CPU or GPU. This leads to significant resource waste on existing heterogeneous computing platforms (where both GPUs and CPUs coexist). For example, during the DC and BC-Back phases, only CPU resources are used, leaving the GPU idle, resulting in GPU resource waste. Another example is the SBR-Back phase, which only uses GPU resources, leaving the CPU idle, resulting in CPU resource waste. Furthermore, the existing BC-Back phase is implemented on the CPU, resulting in insufficient parallelism. Summary of the Invention

[0003] To address the issue of CPU or GPU resource waste in symmetric matrix eigenvalue decomposition and the lack of processor parallelism in the existing convex-chasing backward transform (BC-Back) process on CPUs, we propose a novel and efficient method, system, and device for solving eigenvalues ​​of heterogeneous symmetric matrices based on GPU parallel acceleration.

[0004] This method mainly involves two parts:

[0005] In the first part, a novel and efficient convex block chasing backward transformation method based on GPU parallel acceleration is proposed. It directly uses the Householder vector based on the GPU kernel function to perform efficient parallel calculations. At the same time, it no longer transforms the matrix output by the divide-and-conquer stage, but transforms the matrix output by the SBR-Back stage (the SBR-Back stage no longer transforms the output of BC-Back, but directly transforms the unit matrix I to obtain , Express Perform transpose operation). First, it treats the transformation matrix The data is split into columns and the maximum number of thread blocks that the GPU can create is created. Each thread block processes the columns in parallel. This part maximizes the parallelization of the GPU and maximizes the utilization of the GPU's computing resources.

[0006] Part 2: Since the input of BC-Back in Part 1 of this method no longer depends on the DC stage and it is implemented on the GPU, it can be processed in parallel with the DC stage using a sub-process that is completely executed on the CPU. At the same time, in order to fully utilize resources, the SBR-Back stage is placed before the BC-Back stage, allowing the entire process of these two stages to be parallel to the DC stage. The present invention no longer follows the serial order of SBR -> BC -> DC -> BC-Back -> SBR-Back, but uses the parallelization of the stages consisting of SBR -> BC, DC, and (SBR-Back + BC-Back) in series, and introduces the matrix multiplication GEMM method for calculation, see Figure 1 As shown in the figure below, this method realizes parallelism and fully utilizes the computing resources of the graphics processing unit (GPU) and the central processing unit (CPU), solving the problem of insufficient resource utilization in existing serial processes and significantly improving the processor's processing capabilities in symmetric matrix eigenvalue decomposition.

[0007] The present invention provides a novel and efficient GPU-based parallel acceleration method for heterogeneous symmetric matrix eigenvalue decomposition, which specifically includes the following steps:

[0008] Step 1: Given the original matrix A to be solved, obtain the unit orthogonal matrix, tridiagonal matrix and Householder vector matrix through strip decomposition and tridiagonal decomposition; then solve the tridiagonal matrix to obtain the eigenvalue array of the original matrix and the eigenvector matrix of the tridiagonal matrix;

[0009] Step 2: Using the BC-Back method, perform a Householder transformation on the unit orthogonal matrix using the Householder vectors in the Householder vector matrix to obtain the intermediate process matrix;

[0010] Specifically, a striped array is allocated on the GPU, and the Householder vector u in the Householder vector matrix is ​​striped and stored to improve the efficiency of the GPU in reading the Householder vector. Then, the maximum size of the shared memory that can be used by the thread block on the GPU, the number of streaming multiprocessors, and the number of thread blocks that can be started on the streaming multiprocessor are obtained. Based on the number of streaming multiprocessors and the number of thread blocks that can be started on each streaming multiprocessor, an available thread block is created. Each thread of the available thread block uses the Householder vector to transform the transposed unit orthogonal matrix to obtain an intermediate process matrix.

[0011] Step 3: Use matrix multiplication to multiply the intermediate process matrix and the eigenvector matrix to complete the matrix eigenvalue decomposition.

[0012] The step 2 is specifically as follows:

[0013] Step 2.1: After striping the Householder vector u in the Householder vector matrix, allocate a two-dimensional shared memory array on the GPU to store the Householder vector.

[0014] Step 2.2: Copy the Householder vector from the Householder vector matrix to the two-dimensional shared memory array;

[0015] Step 2.3: The transposed unit orthogonal matrix is ​​used as the matrix to be transformed, and each thread creates an array r q Store the columns of the matrix to be transformed assigned to this thread, for array r q Perform Householder transformation;

[0016] Step 2.4: Each thread traverses all Householder vectors and completes the transformation operation.

[0017] Furthermore, a possible implementation of step 2 is as follows:

[0018] Assume that the size of the original matrix A is n*n, the length of the Householder vector u is b, and the number of thread warps in each thread block on the GPU is k. Solve the intermediate process matrix :

[0019] Step 2.1: Obtain the number of Householder vectors u; obtain the total number of bump chases, denoted as sN, where sN = (n + b - 3) / b; obtain the number of Householder vectors u generated by the last bump chase, denoted as tUN, calculated using the formula tUN = n - (sN - 1) * b - 2;

[0020] Step 2.2: Allocate and store the striped array bandUArrey on the GPU, and store the Householder vector u stored in the Householder vector matrix into the array bandUArrey;

[0021] Step 2.3: Create a thread block; get the maximum size of shared memory that can be used by the thread block maxSharedMenPerBlock, the number of stream multiprocessors on the GPU, and the number of thread blocks that can be started on the stream multiprocessor; set the size of shared memory available in the thread block to for:

[0022]

[0023] Among them, Sizeof() is a function in C language that finds the memory length occupied by the data type datatype; let the maximum number of vectors u that can be loaded by the thread block at one time maxSharedUNum be:

[0024] maxSharedUNum= sharedMemSize / (sizeof(datatype)*b)

[0025] The number of threads in a thread block dimblock is (32, k, 1), where 32 is the number of threads in a warp in an NVIDIA GPU.

[0026] Create blockNum thread blocks based on the number of stream multiprocessors and the number of thread blocks that can be started on each stream multiprocessor;

[0027] Step 2.4: Transformation begins; each thread transposes the identity orthogonal matrix The perBlockN column of the array bandUArrey is traversed for transformation, where perBlockN=(n+blockNum-1) / blockNum, where blockNum is the number of thread blocks created and perBlockN is the index of the column to be transformed;

[0028] Step 2.5: Transformation is completed; wait for all thread blocks to complete the transformation operation and obtain the intermediate process matrix.

[0029] Furthermore, the step 2.4 is specifically as follows:

[0030] Step 2.4.1: Allocate a two-dimensional shared memory array sU[maxSharedUNum][b] on the GPU to store the vector u;

[0031] Step 2.4.2: Get the total number of vectors u in this bump chase, procTotalUNum = i*b + tUN; set the number of vectors u that can be loaded into the two-dimensional shared memory array sU, procUNum = min(maxSharedUNum, procTotalUNum-j), and load procUNum u from the striped array bandUArrey[indexU] into the array sU;

[0032] Step 2.4.3: Each thread block loads the corresponding Column to array r q In, r q It is a memory array created on the GPU to store the Column; for each thread block r q Perform Householder transformation:

[0033] r q =(I-sU[t]*sU[t] T )* r q ;

[0034] Where I is the identity matrix.

[0035] The present invention also provides an application of a heterogeneous symmetric matrix eigenvalue decomposition method based on GPU parallel acceleration in image compression, which specifically includes the following steps:

[0036] Step 1: Read the original image to be compressed, call the conversion function to convert it into an RGB image, and then store the RGB red, green, and blue channel data in three two-dimensional array matrices of red, green, and blue respectively. In the red, green and blue arrays Pass it to the GPU and call the function on these three array matrices , solve the three matrix multiplications and get three symmetric matrices for red, green, and blue:

[0037]

[0038]

[0039]

[0040] Step 2: For the three symmetric matrices of red, green and blue Call the symmetric matrix eigenvalue decomposition method proposed in this invention to obtain The corresponding eigenvalue array 、 、 and the eigenvector matrix 、 、 ; The details are as follows:

[0041] Step 2.1: Take the red symmetric matrix For example, through strip decomposition and triangular decomposition, we can get the matrix composed of Householder vectors in the SBR process: , a tridiagonal matrix and a matrix consisting of Householder vectors in the BC process ; Solve the tridiagonal matrix and get the red symmetric matrix The eigenvalue array of The eigenvector matrix corresponding to the tridiagonal matrix .

[0042] Step 2.2: For the identity matrix The matrix composed of Householder vectors in the SBR process The Householder vector in is transformed by Householder to obtain the unit orthogonal matrix ;

[0043] Step 2.3: Using the BC-Back method, the unit orthogonal matrix Use the matrix composed of Householder vectors in the BC process The Householder vector in is transformed by Householder to obtain the intermediate process matrix ;

[0044] Step 2.4: Use matrix multiplication to convert the intermediate process matrix and the eigenvector matrix Multiply to get the eigenvector matrix , complete the matrix eigenvalue decomposition;

[0045] Similarly, the green and blue symmetric matrices After steps 2.1-2.4, the symmetric matrices of green and blue can be obtained The eigenvalue array of 、 and the eigenvector matrix 、 ;

[0046] Step 3: Get the eigenvalue array , , The number of elements with the first set percentage as the eigenvalue , , ; Extract the eigenvector matrix , , The first set percentage of the column as the final eigenvector matrix , , ;

[0047] Step 4: Red, green, and blue arrays , , Call the square root function , find the square root value , , ;

[0048] Step 5: Call the function , calculate the matrix multiplication respectively, and get three matrix arrays:

[0049]

[0050]

[0051]

[0052] Step 6: Final eigenvector matrix of red, green and blue channels 、 、 , matrix array 、 、 and eigenvalues 、 、 Copy it back to the CPU memory and store it on disk to complete the image compression.

[0053] The present invention also provides a heterogeneous symmetric matrix eigenvalue decomposition system based on GPU parallel acceleration to implement the above method, the system comprising:

[0054] The original matrix solving unit obtains the unit orthogonal matrix, tridiagonal matrix and Householder vector matrix through strip decomposition and tridiagonal decomposition; solves the tridiagonal matrix to obtain the eigenvector matrix corresponding to the tridiagonal matrix;

[0055] GPU bump chase backpropagation unit, transpose of the unit orthogonal matrix Perform Householder transformation; by allocating memory arrays on the GPU, create available thread blocks based on the number of stream multiprocessors on the GPU and the number of thread blocks that can be started on each stream multiprocessor, and each thread in the available thread block transforms the transposed unit orthogonal matrix using the Householder vector in the Householder vector matrix to obtain an intermediate process matrix;

[0056] The GEMM unit is used to multiply the intermediate process matrix and the eigenvector matrix to complete the matrix eigenvalue decomposition.

[0057] The present invention also provides a heterogeneous symmetric matrix eigenvalue decomposition device based on GPU parallel acceleration, including a memory, a processor and a computer program stored in the memory, the processor including a central processing unit CPU and a graphics processing unit GPU, and the processor executes the computer program to implement the steps of the above method.

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

[0059] The present invention modifies the BC-Back action matrix in the traditional symmetric matrix eigenvalue decomposition method to be an SBR-Back output matrix (the traditional action matrix is ​​a DC output matrix), breaking the design that the five sub-processes of the symmetric matrix eigenvalue decomposition in the prior art can only be strictly executed in series, and proposes a novel CPU / GPU parallel processing method for solving the symmetric matrix eigenvalues, solving the problem of idle CPU and GPU computing resources when the five sub-processes are serial, and greatly improving the computing efficiency. At the same time, the present invention is aimed at the BC-Back stage and proposes a novel BC-Back method based on GPU implementation, which transfers the BC-Back process from a CPU with limited parallel capability to a GPU that can be large-scale parallel, greatly improving the parallel efficiency of the BC-Back stage. The present invention makes full use of the characteristics of heterogeneous architecture and GPU parallelism, realizes the parallelism between the symmetric matrix eigenvalue solving sub-processes and the internal parallelism of BC-Back, makes full use of CPU and GPU resources, and greatly improves the computing performance of the computer for the symmetric matrix eigenvalue solving process.

[0060] In an environment consisting of an NVIDIA H100 GPU and an Intel(R) Xeon(R) Platinum 8462Y CPU, for a matrix size of 32768*32768, this method achieves a 2.82-fold performance improvement in stage 4 compared to the existing symmetric matrix eigenvalue decomposition method MAGAM. The overall symmetric matrix eigenvalue decomposition achieves a 2.51-fold and 2.38-fold performance improvement compared to MAGMA and CUDA cuSolver, respectively.

[0061] The method provided by the present invention can also be applied to the field of image compression, and while ensuring the image compression rate, it improves the processing speed during image compression and ensures a smaller degree of distortion during image compression. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 2 is a flow chart comparing the method of the present invention and the traditional method. DETAILED DESCRIPTION

[0063] In order to better understand the purpose, structure and function of the present invention, the following is a further detailed description of a method for solving eigenvalues ​​of heterogeneous symmetric matrices based on GPU parallel acceleration of the present invention in conjunction with the accompanying drawings.

[0064] Example 1

[0065] A method for eigenvalue decomposition of heterogeneous symmetric matrices based on GPU parallel acceleration, the specific steps are as follows:

[0066] Step 1: Given the original matrix A to be solved, obtain the unit orthogonal matrix, tridiagonal matrix and Householder vector matrix through strip decomposition and tridiagonal decomposition; then solve the tridiagonal matrix to obtain the eigenvalue array of the original matrix and the eigenvector matrix of the tridiagonal matrix;

[0067] Step 2: Using the BC-Back method, perform a Householder transformation on the unit orthogonal matrix using the Householder vectors in the Householder vector matrix to obtain the intermediate process matrix;

[0068] Step 2.1: Allocate a striped array on the GPU and store the Householder vector u in the Householder vector matrix in a striped manner to improve the GPU's reading efficiency of the Householder vector.

[0069] Step 2.2: Obtain the maximum size of the shared memory available to thread blocks on the GPU, the number of streaming multiprocessors, and the number of thread blocks that can be launched on each streaming multiprocessor; create an available thread block based on the number of streaming multiprocessors and the number of thread blocks that can be launched on each streaming multiprocessor. Each thread in the available thread block transforms the transposed unit orthogonal matrix using a Householder vector to obtain an intermediate process matrix;

[0070] Step 3: Use matrix multiplication to multiply the intermediate process matrix and the eigenvector matrix to complete the matrix eigenvalue decomposition.

[0071] The step 2.2 is as follows:

[0072] Step 2.2.1: Allocate a two-dimensional shared memory array on the GPU to store the Householder vector;

[0073] Step 2.2.2: Copy the Householder vector from the Householder vector matrix to the two-dimensional shared memory array;

[0074] Step 2.2.3: The transposed unit orthogonal matrix is ​​used as the matrix to be transformed, and each thread creates an array rq Store the columns of the matrix to be transformed assigned to this thread, for array r q Perform Householder transformation;

[0075] Step 2.2.4: Each thread traverses all Householder vectors and completes the transformation operation.

[0076] Example 2

[0077] In step 2 of Example 1, the intermediate process matrix is ​​obtained using the BC-Back method. Another possible implementation method is as follows:

[0078] Assume that the size of the original matrix A is n*n, the length of the Householder vector u is b, and the number of thread warps in each thread block on the GPU is k. Solve the intermediate process matrix :

[0079] Step 2.1: Get the total number of bump chases, denoted as sN, where sN = (n + b - 3) / b. Get the number of Householder vectors u generated by the last bump chase, denoted as tUN, calculated using the formula tUN = n - (sN - 1) * b - 2.

[0080] Step 2.2: Allocate a striped array bandUArrey on the GPU and store the Householder vector u stored in the Householder vector matrix into the array bandUArrey;

[0081] Step 2.3: Get the maximum size of shared memory that can be used by the thread block, maxSharedMemPerBlock, the number of streaming multiprocessors on the GPU, and the number of thread blocks that can be started on the streaming multiprocessors; let the size of shared memory available in the thread block be for:

[0082]

[0083] Among them, Sizeof() is a function in C language that finds the memory length occupied by the data type datatype; let the maximum number of vectors u that can be loaded by the thread block at one time maxSharedUNum be:

[0084] maxSharedUNum= sharedMemSize / (sizeof(datatype)*b)

[0085] The number of threads in a thread block dimblock is (32, k, 1), where 32 is the number of threads in a warp in an NVIDIA GPU.

[0086] Create blockNum thread blocks based on the number of stream multiprocessors and the number of thread blocks that can be started on each stream multiprocessor;

[0087] Step 2.4: Each thread transposes the identity orthogonal matrix The perBlockN column of the array bandUArrey is traversed for transformation, where perBlockN=(n+blockNum-1) / blockNum;

[0088] Step 2.5: Wait for all thread blocks to complete the transformation operation and obtain the intermediate process matrix.

[0089] Furthermore, the step 2.4 is specifically as follows:

[0090] Step 2.4.1: Allocate a two-dimensional shared memory array sU[maxSharedUNum][b] on the GPU to store the vector u;

[0091] Step 2.4.2: Get the total number of vectors u in this bump chase, procTotalUNum = i*b + tUN; set the number of vectors u that can be loaded into the two-dimensional shared memory array sU, procUNum = min(maxSharedUNum, procTotalUNum-j), and load procUNum u from the striped array bandUArrey[indexU] into the array sU;

[0092] Step 2.4.3: Each thread block loads the corresponding Column to array r q In, r q It is a memory array created on the GPU to store the Column; for each thread block r q Perform Householder transformation:

[0093] r q =(I-sU[t]*sU[t] T )* r q ;

[0094] Where I is the identity matrix.

[0095] Example 3

[0096] Taking NVIDIA's H100 GPU as an example, assuming that the size of the original matrix A for eigenvalue decomposition is n*n, the matrix dimension n=32768, the length of the Householder vector u b=32, the number of thread warps in each thread block k=32, and the data type to be processed dataType=double, the specific steps of the heterogeneous symmetric matrix eigenvalue decomposition method based on GPU parallel acceleration provided in this embodiment are as follows:

[0097] Step 1: Complete the striped decomposition of the original matrix A to be solved, which is expressed in the following formula:

[0098] A=Q1*B*Q1 T

[0099] Where A is the original matrix, B is the striped matrix obtained after striped decomposition, and Q1 is the unit orthogonal matrix. However, in general, Q1 is not the unit orthogonal matrix. In order to reduce the amount of calculation, the WY representation is generally used to store two process matrices W and Y.

[0100] Step 2: Decompose the striped matrix B into tridiagonal form. The formula is as follows:

[0101] B=Q2*T*Q2 T

[0102] That is, B=U*T*U T

[0103] Among them, T is the tridiagonal matrix obtained by tridiagonal decomposition. Similarly, to reduce the amount of calculation, the matrix Q2 does not store the unit orthogonal matrix, but the Householder vector u generated by the tridiagonal decomposition. For the purpose of unified description later, the matrix U is equivalent to the matrix Q2.

[0104] Step 3: Use the Divide and Conquer (DC) algorithm in the MKL library to solve the tridiagonal matrix T. The specific formula is as follows:

[0105] T=Q T *Λ* Q T T

[0106] where Q T is the eigenvector matrix corresponding to the tridiagonal matrix T, and the matrix Λ is a diagonal matrix whose diagonal elements are the eigenvalues ​​of the original matrix A.

[0107] At the same time, the existing SBR-Back method on the GPU (such as the SBR-Back method in MAGMA) is used in parallel to transform the intermediate process matrices W and Y to obtain the unit orthogonal matrix Q1, and the transpose operation is performed on it to obtain Q1 T Then use the BC-Back method to find the intermediate process matrix Q4. The BC-Back method formula is simplified as follows:

[0108] Q4=Q2 T *Q1 T

[0109] Step 4: Finally, use matrix multiplication to find the unit orthogonal matrix Q = Q4 T *Q T T , complete the matrix eigenvalue decomposition, and the matrix Q is the eigenvector matrix of the original matrix A.

[0110] Furthermore, the BC-Back method obtains the intermediate process matrix Q4 through the following steps:

[0111] Step 3.1: Calculate the total number of Householder vectors u in the matrix U and denote it as tN. Let tN=0. The total number of bump chases is sN, sN=(n+b-3) / b, which is (32768+32-3) / 32. The vector generated by the last bump chase is The number of tUN=n-(sN-1)*b-2, that is, 32768*(sN-1)*32-2; the first iteration variable i=0, jump to step 3.1.1;

[0112] Step 3.1.1: Determine whether i is less than sN. If not, jump to step 3.2; otherwise, tN=tN+32*i+tUN, set i=i+1, and jump to step 3.1.1;

[0113] Step 3.2: Allocate a two-dimensional array bandUArrey[tN]

[32] on the GPU for striped storage of vector u in matrix U; let i = 0, index indexU of vector u = 0; jump to step 3.2.1;

[0114] Step 3.2.1: Determine whether i is less than sN. If not, jump to step 3.3; if less, jump to step 3.2.2.

[0115] Step 3.2.2: Set the second iteration variable j = 0; the starting position of the vector u to be loaded (stored) is base = (sN-(i+1))*32+1;

[0116] Let the total number of vectors u to be loaded in this process be h=32*i+tUN; jump to step 3.2.3;

[0117] Step 3.2.3: Determine whether j is less than the total number h of vectors u. If not, set i++ and jump to step 3.2.1. If so, copy u to the striped array bandUArrey. The pseudo code is as follows:

[0118] memCpy(U+base+j*n,bandUArrey[indexU]);

[0119] Then let indexU++, j++; jump to step 3.2.3;

[0120] Step 3.3: Calculate the number of thread blocks that can be created. Assume that the element data type of matrix U is double. Use the GPU library function to obtain the maximum size of shared memory that can be used by thread blocks, maxSharedMemPerBlock, the number of SMs (streaming multiprocessors) on the GPU, SMCount, and the number of thread blocks that can be launched on the SM, numBlockPerSm. Let the size of shared memory available in the thread block be:

[0121]

[0122] Among them, Sizeof() is a function in C language to find the memory length occupied by double type data; let the maximum number of vectors u that can be loaded by the thread block at one time, maxSharedUNum, be:

[0123] maxSharedUNum= sharedMemSize / (sizeof(double)*32)

[0124] The number of threads in the thread block dimblock is (32,32,1);

[0125] Create blockNum thread blocks based on the number of SMs and the number of thread blocks that can be started on each SM, blockNum=SMCount*numBlocksPerSm; each thread is responsible for Q1 T The perBlockN column of the matrix is ​​transformed using the vector u in the matrix U, perBlockN=(32768+blockNum-1) / blockNum; the specific process is as follows:

[0126] Step 3.3.1: Get the thread block index blockIndex of each thread block; the starting column number for each thread to process startCol=blockIndex*perBlockN; determine whether startCol is less than n, if not, jump to step 3.4; if so, set perBlockN=min(perBlockN,n-startCol), i=0, indexU=0, allocate a two-dimensional shared memory array sU[maxSharedUNum]

[32] on the GPU to store the vector u, and jump to step 3.3.2;

[0127] Step 3.3.2: Determine whether i is less than sN. If not, jump to step 3.4. If so, set j = 0. The total number of vectors u in this bump chase is procTotalUNum = i * 32 + tUN. Jump to step 3.3.3.

[0128] Step 3.3.3: Determine whether j is less than procTotalUNum. If not, set i = i + 1 and jump to step 3.3.2. Otherwise, set procUNum = min(maxSharedUNum, procTotalUNum - j) for the number of vectors u that can be loaded into the two-dimensional array sU. Load procUNum u from the array bandUArrey[indexU] into the array sU. Set h = 0 and jump to step 3.3.4.

[0129] Step 3.3.4: Determine whether h is less than perBlockN. If not, set j = j + procUNum, indexU = indexU + procUNum, and jump to step 3.3.3; otherwise, set Q1 that can be processed by one thread block at a time. T The number of columns strideN=min(k, perBlockN-h), load strideN columns Q1 T Processing; the specific steps are as follows:

[0130] Step 3.3.4.1: Each thread warp loads the corresponding Q1 T Column to array r q

[32] In q

[32] is the array created on the GPU to store the Q1 to be processed T Column, set the third iteration variable t=0, jump to step 3.3.4.2;

[0131] Step 3.3.4.2: Determine whether t is less than procUNum. If not, set r q

[32] Save back the corresponding Q1 TIn the example, let h = h + strideN and jump to step 3.3.4; otherwise, for each thread bundle r q

[32] Perform Householder transformation:

[0132] r q =(I-sU[t]*sU[t] T )* r q

[0133] Where I is the identity matrix, set t = t + 1, and jump to step 3.3.4.2;

[0134] Step 3.4: Wait for all thread blocks to complete the operation and get the latest Q1 T This is the matrix Q4 to be determined.

[0135] A novel and efficient method for solving the eigenvalues ​​of heterogeneous symmetric matrices based on GPU parallel acceleration is provided in the present embodiment. First, a novel BC-Back method based on GPU acceleration is used in the BC-Back stage. It is all implemented on the GPU, making full use of the performance of the existing GPU architecture, and can bring higher computational parallelism than the prior art. In addition, its input no longer depends on the output of the DC stage, but uses the output of the SBR-Back stage as input, so the BC-Back stage can be carried out in parallel with the DC stage. According to the improvement in the BC-Back stage, the present invention parallelizes the DC stage of the whole eigenvalue solving process and the entirety of the (SBR-Back+BC-Back) serial composition. Since the DC stage is executed on the CPU, the SBR-Back stage+BC-Back stage is executed on the GPU. The two parallel processes fully utilize the performance of the existing CPU / GPU heterogeneous architecture, solve the problem that the GPU / CPU idleness existing in the prior art causes waste of resources, and make full use of current GPU and CPU resources.

[0136] The proposed method, in an environment consisting of an NVIDIA H100 GPU and an Intel(R) Xeon(R) Platinum 8462Y CPU, achieves a 2.82-fold performance improvement in the BC-Back stage compared to the existing MAGAM method when the matrix size is 32768*32768, and the entire symmetric matrix eigenvalue decomposition achieves a 2.51-fold and 2.38-fold performance improvement compared to MAGMA and cuSolver, respectively. When the matrix size is 45056*45056, stage 4 achieves a 2.75-fold performance improvement compared to MAGAM, and the entire symmetric matrix eigenvalue decomposition (EVD) achieves a 2.54-fold and 2.59-fold performance improvement compared to MAGMA and cuSolver, respectively, as shown in Table 1.

[0137] Table 1: Comparison of eigenvalue decomposition of symmetric matrices using different methods

[0138] Method Name Matrix size BC-Back phase (unit: s) The entire EVD process (unit: s) MAGMA 32768 22.13 42.62s cuSolver 32768 \ 40.5 This method 32768 7.86 17.00s MAGMA 45056 55.55 99.3 cuSolver 45056 \ 101.4 This method 45056 20.21 39.08

[0139] Example 4

[0140] This embodiment further provides a heterogeneous symmetric matrix eigenvalue decomposition system based on GPU parallel acceleration to implement the above method, the system comprising:

[0141] The original matrix solving unit obtains the unit orthogonal matrix, tridiagonal matrix and Householder vector matrix through strip decomposition and tridiagonal decomposition; solves the tridiagonal matrix to obtain the eigenvector matrix corresponding to the tridiagonal matrix;

[0142] GPU bump chase backpropagation unit, transpose of the unit orthogonal matrix Perform Householder transformation; by allocating memory arrays on the GPU, create available thread blocks based on the number of stream multiprocessors on the GPU and the number of thread blocks that can be started on each stream multiprocessor, and each thread in the available thread block transforms the transposed unit orthogonal matrix using the Householder vector in the Householder vector matrix to obtain an intermediate process matrix;

[0143] The GEMM unit is used to multiply the intermediate process matrix and the eigenvector matrix to complete the matrix eigenvalue decomposition.

[0144] Example 5

[0145] This embodiment also provides a heterogeneous symmetric matrix eigenvalue decomposition device based on GPU parallel acceleration, including a memory, a processor and a computer program stored in the memory, the processor including a central processing unit CPU and a graphics processing unit GPU, and the processor executes the computer program to implement the steps of the above method.

[0146] Example 6

[0147] This embodiment provides an application of a heterogeneous symmetric matrix eigenvalue decomposition method based on GPU parallel acceleration in image compression, which specifically includes the following steps:

[0148] Step 1: Read the original image to be compressed, call the conversion function to convert it into an RGB image, and then store the RGB three-channel data in a two-dimensional array matrix. In this example, we will use these three arrays Pass it to the GPU and call the function on these three array matrices , solve the three matrix multiplications and get three symmetric matrices:

[0149]

[0150]

[0151]

[0152] Step 2: Symmetrically transform the three matrices Call the symmetric matrix eigenvalue decomposition method proposed in this invention to obtain The corresponding eigenvalue array 、 、 and the eigenvector matrix 、 、 ; The details are as follows:

[0153] Step 2.1: Symmetric matrix of channels For example, through strip decomposition and triangular decomposition, we can get the matrix composed of Householder vectors in the SBR process: , a tridiagonal matrix and a matrix consisting of Householder vectors in the BC process ; Solve the tridiagonal matrix and get Channel symmetric matrix The eigenvalue array of The eigenvector matrix corresponding to the tridiagonal matrix .

[0154] Step 2.2: For the identity matrix The matrix composed of Householder vectors in the SBR process The Householder vector in is transformed by Householder to obtain the unit orthogonal matrix ;

[0155] Step 2.3: Using the BC-Back method, the unit orthogonal matrix Use the matrix composed of Householder vectors in the BC process The Householder vector in is transformed by Householder to obtain the intermediate process matrix ;

[0156] Step 2.4: Use matrix multiplication to convert the intermediate process matrix and the eigenvector matrix Multiply to get the eigenvector matrix , complete the matrix eigenvalue decomposition;

[0157] Similarly, the symmetric matrix After steps 2.1-2.4, the symmetric matrix can be obtained The eigenvalue array of 、 and the eigenvector matrix 、 ;

[0158] Step 3: Get the eigenvalue array , , The number of elements with the first set percentage as the eigenvalue , , ; Extract the eigenvector matrix , , The first set percentage of the column as the final eigenvector matrix , , ;

[0159] Step 4: Array , , Call the square root function , find the square root value , , ;

[0160] Step 5: Call the function , calculate the matrix multiplication respectively, and get three matrix arrays:

[0161]

[0162]

[0163]

[0164] Step 6: Final eigenvector matrix 、 、 , matrix array 、 、 and eigenvalues 、 、 Copy it back to the CPU memory and store it on disk to complete the image compression.

[0165] By applying this method to image compression, the image compression ratio is approximately 6.4x, and the image distortion is less than 5%. By applying the proposed symmetric matrix eigenvalue decomposition method based on a heterogeneous GPU and CPU parallel processing architecture, this novel image compression method can increase the compression speed of a single 720p image by approximately 1.4x compared to the AVIF compression method.

[0166] It will be understood that the present invention is described by way of some embodiments, and it will be appreciated by 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 to the specific embodiments disclosed herein, and all embodiments falling within the scope of the claims of this application are intended to be protected by the present invention.

Claims

1. A method for eigenvalue decomposition of heterogeneous symmetric matrices based on GPU parallel acceleration, characterized in that: The method runs in parallel on a graphics processing unit (GPU) and a central processing unit (CPU), and specifically includes the following steps: Step 1: Given the original matrix A to be solved, obtain the unit orthogonal matrix, tridiagonal matrix, and Householder vector matrix through strip decomposition and tridiagonal decomposition; then solve the tridiagonal matrix to obtain the eigenvalue array of the original matrix and the eigenvector matrix of the tridiagonal matrix; Step 2: Using the convex block chasing backward transformation method, the unit orthogonal matrix is ​​transformed using the Householder vector in the Householder vector matrix to obtain the intermediate process matrix; Specifically, first, a striped array is allocated on the GPU, and the Householder vector u in the Householder vector matrix is ​​striped and stored; then, the maximum size of the shared memory available to the thread block on the GPU, the number of streaming multiprocessors, and the number of thread blocks that can be started on the streaming multiprocessors are obtained; based on the number of streaming multiprocessors and the number of thread blocks that can be started on each streaming multiprocessor, an available thread block is created, and each thread in the available thread block uses the Householder vector to transform the transposed unit orthogonal matrix to obtain an intermediate process matrix; assuming that the size of the original matrix A is n*n, the length of the Householder vector u is b, and the number of thread warps in each thread block on the GPU is k, the solution process of the intermediate process matrix is ​​as follows: Get the number of Householder vectors u; get the total number of bump chases, denoted as sN, where sN = (n + b - 3) / b; get the number of Householder vectors u generated by the last bump chase, denoted as tUN, calculated using the formula tUN = n - (sN - 1) * b - 2; Allocation and storage of striped arrays; allocating a striped array bandUArrey on the GPU, and storing the Householder vector u stored in the Householder vector matrix into the striped array bandUArrey; Create a thread block; get the maximum size of shared memory that can be used by the thread block, maxSharedMemPerBlock, the number of streaming multiprocessors on the GPU, and the number of thread blocks that can be started on the streaming multiprocessors; let the size of shared memory available in the thread block, sharedMemSize, be: Among them, Sizeof() is a function in C language that finds the memory length occupied by the data type datatype; let the maximum number of vectors u that can be loaded by the thread block at one time maxSharedUNum be: maxSharedUNum=sharedMemSize / (sizeof(datatype)*b) The number of threads in the thread block dimblock is (32, k, 1), where 32 is the number of threads in the thread warp in the GPU; blockNum thread blocks are created based on the number of streaming multiprocessors and the number of thread blocks that can be started on each streaming multiprocessor; the transformation begins; each thread transposes the unit orthogonal matrix The perBlockN column of the striped array bandUArrey is traversed for transformation, where perBlockN = (n + blockNum - 1) / blockNum, where blockNum is the number of thread blocks created and perBlockN is the index of the column to be transformed; The transformation is completed; wait for all thread blocks to complete the transformation operation and obtain the intermediate process matrix; Step 3: Use matrix multiplication to multiply the intermediate process matrix and the eigenvector matrix to complete the matrix eigenvalue decomposition.

2. The method for eigenvalue decomposition of heterogeneous symmetric matrices based on GPU parallel acceleration according to claim 1, characterized in that: The step 2 is specifically as follows: Step 2.1: Stripe the Householder vector u in the Householder vector matrix and allocate a two-dimensional shared memory array on the GPU to store the Householder vector. Step 2.2: Copy the Householder vectors from the Householder vector matrix to a two-dimensional shared memory array. Step 2.3: The transposed unit orthogonal matrix is ​​used as the matrix to be transformed, and each thread creates an array r q Store the columns of the matrix to be transformed assigned to this thread, for array r q Perform Householder transformation; Step 2.4: Each thread traverses all Householder vectors and completes the transformation operation.

3. The method for eigenvalue decomposition of heterogeneous symmetric matrices based on GPU parallel acceleration according to claim 2, characterized in that: The transformation starts specifically including the following sub-steps: First, allocate a two-dimensional shared memory array sU[maxSharedUNum][b] on the GPU to store the vector u; Then, obtain the total number of vectors u in this bump chase: procTotalUNum = i*b + tUN; set the number of vectors u that can be loaded into the two-dimensional shared memory array sU: procUNum = min(maxSharedUNum, procTotalUNum - j); load procUNum u from the striped array bandUArrey[indexU] into the two-dimensional shared memory array sU; Finally, each thread block loads the corresponding Column to array r q In, r q It is a memory array created on the GPU to store the Column; for each thread block r q Perform a Householder transform: r q =(I-sU[t]*sU[t] T )*r q ; Where I is the identity matrix.

4. The method for eigenvalue decomposition of heterogeneous symmetric matrices based on GPU parallel acceleration according to any one of claims 1 to 3, characterized in that: Applied to image compression, specifically including the following steps: Step S1: Read the original image to be compressed, call the conversion function to convert it into an RGB image, and then store the RGB red, green and blue channel data in three two-dimensional array matrices R, G, and B respectively; transfer the three red, green and blue arrays R, G, and B to the GPU, call the GEMM function on these three array matrices, solve the three matrix multiplications, and obtain three symmetric matrices: A R =RR T ; A G =GG T ; A B =BB T ; Step S2: For the three symmetric matrices A for red, green and blue R , A G , A B Call a heterogeneous symmetric matrix eigenvalue decomposition method based on GPU parallel acceleration to obtain A R 、A G 、A B The corresponding red, green, and blue eigenvalue arrays ∑ R ,∑ G ,∑ B and the red, green, and blue eigenvector matrices Q R , Q G , Q B ; The details are as follows: Step S2.1: Take the red symmetric matrix A R For example, through strip decomposition and triangular decomposition, we can get the matrix U composed of the Householder vectors in the SBR process. s , a tridiagonal matrix and a matrix U consisting of Householder vectors in the BC process b ; Solve the tridiagonal matrix and get the red symmetric matrix A R The eigenvalue array ∑ R The eigenvector matrix Q corresponding to the tridiagonal matrix d ; Step S2.2: Apply the Householder vector matrix U to the identity matrix I s Perform Householder transformation on the Householder vector in to obtain the unit orthogonal matrix Q s ; Step S2.3: Using the convex block chasing backward transformation method, the unit orthogonal matrix Q s The matrix U composed of Householder vectors in the BC process is used b Perform Householder transformation on the Householder vector in to obtain the intermediate process matrix Q1; Step S2.4: Use matrix multiplication to combine the intermediate process matrix Q1 and the eigenvector matrix Q d Multiply to get the eigenvector matrix Q R , complete the matrix eigenvalue decomposition; Similarly, the green and blue symmetric matrices A G 、A B After steps 2.1-2.4, we can get the green and blue symmetric matrix A G 、A B The eigenvalue array ∑ G ,∑ B and the eigenvector matrix Q G , Q B ; Step S3: Take out the red, green and blue eigenvalue arrays ∑ R ,∑ G ,∑ B The number of elements in the first set percentage is used as the array eigenvalue σ R , σ G , σ B ; Take out the red, green and blue eigenvector matrix Q R , Q G , Q B The first set percentage of columns is used as the final eigenvector matrix U R , U G , U B ; Step S4: Calculate the red, green and blue array eigenvalues ​​σ R , σ G , σ B Call the square root calculation function sqrt to find its square root value σ′ R =sqrt(σ R ),σ′ G =sqrt(σ G ),σ′ B =sqrt(σ B ); Step S5: Call the GEMM function to calculate matrix multiplication and obtain three matrix arrays: red, green, and blue: Step S6: The final eigenvector matrix U of the red, green and blue channels R 、U G 、U B , matrix array V R 、V B 、V G and array eigenvalue σ R , σ G , σ B Copy it back to the CPU memory and store it on disk to complete the image compression.

5. A heterogeneous symmetric matrix eigenvalue decomposition system based on GPU parallel acceleration, characterized in that: The system is used to implement the method according to any one of claims 1 to 3, and the system comprises: The original matrix solving unit obtains the unit orthogonal matrix, tridiagonal matrix and Householder vector matrix through striping decomposition and tridiagonal decomposition; solves the tridiagonal matrix to obtain the eigenvector matrix corresponding to the tridiagonal matrix; GPU bump chase back propagation unit, the matrix Perform a Householder transform; by allocating a memory array on the GPU, create an available thread block based on the number of stream multiprocessors on the GPU and the number of thread blocks that can be started on each stream multiprocessor, and each thread in the available thread block transforms the transposed unit orthogonal matrix using the Householder vector in the Householder vector matrix to obtain an intermediate process matrix; The GEMM unit is used to multiply the intermediate process matrix and the eigenvector matrix to complete the matrix eigenvalue decomposition.

6. A device for eigenvalue decomposition of heterogeneous symmetric matrices based on GPU parallel acceleration, comprising a memory, a processor, and a computer program stored in the memory, characterized in that: The processor includes a central processing unit (CPU) and a graphics processing unit (GPU), and the processor executes the computer program to implement the steps of the method according to any one of claims 1 to 3.

Images