A Normalization Method for Eigenvalue Problems of Generalized Hermitian Matrix

Abstract

The present invention relates to a method for normalizing generalized Hermitian matrix eigenvalue problems. This method utilizes parallel heterogeneous computing to combine the Cholesky decomposition of matrix B with the normalization of matrix A. Two streams on the device control the Cholesky decomposition of matrix B and the normalization of matrix A. This allows the device to perform normalization calculations on matrix A while the host is performing work related to matrix B. Furthermore, while the host is performing normalization work on matrix A, the device can also perform work related to matrix B, thus achieving asynchronous operation between the host and device. Specifically, four forms of host-device overlap are achievable: overlap between host-side computation and data transmission between the host and device; overlap between host-side communication and data transmission between the host and device; overlap between device-side computation and data transmission between the host and device; and overlap between host-side communication and device-side computation.

Description

Technical Field

[0001] The present invention relates to the field of high performance computing, and in particular to a method for standardizing generalized Hermitian matrix eigenvalue problems. Background Art

[0002] Solving the eigenvalue and eigenvector problems of large dense matrices is a core computational problem in many areas of computational science and numerical simulation. For example, there are numerous numerical solutions to generalized dense Hermitian eigenvalue problems in electronic structure, material design, simulation, and quantum chemistry calculations. The generalized dense Hermitian eigenvalue problem described in this invention can be described as follows:

[0003] AX=BXΛ (1)

[0004] Wherein, matrices A and B are complex Hermitian dense matrices, matrix B is a positive definite matrix, Λ is a diagonal eigenvalue matrix, and X is an eigenvector matrix. In the prior art, this problem is solved by first converting formula (1) into a standard dense eigenvalue problem and then solving the converted standard eigenvalue problem. The conversion of a generalized Hermitian eigenvalue problem into a standard eigenvalue problem is called a standardization method, which is the key method for solving the generalized eigenvalue problem.

[0005] The current normalization method is based on a two-step technique. In the first step, the matrix B is subjected to Cholesky decomposition: B = LL H , where L is a lower triangular matrix, L H is the conjugate transposed matrix of L. The second step is to calculate the matrix corresponding to the standard eigenvalue problem At this point, the original problem is converted into a standard dense feature problem: Where Y = L H X.

[0006] Current standardized approaches, when used in multi-GPU accelerated computing environments, suffer from the following issues: There's a significant amount of repeated communication between the first and second steps; host-device data transfer and inter-host communication can't work together; inter-host communication and device data processing can't be synchronized; and host- and device-level computations can't be performed simultaneously. (Here, the host represents the CPU, and the device represents the GPU.) Summary of the Invention

[0007] The purpose of the present invention is to provide a generalized Hermitian matrix characteristic problem standardization method, aiming to solve the problems existing in the prior art and improve the efficiency of standardization.

[0008] To solve the above problems, the present invention provides a generalized Hermitian matrix eigenvalue problem normalization method, which is applied to a computing system, wherein the computing system includes a host end and a device end, and the device end includes a first stream and a second stream. The method includes:

[0009] The first stream receives a matrix B, and the second stream receives a matrix A, wherein the matrix A is a Hermitian matrix and the matrix B is a positive definite Hermitian matrix;

[0010] The matrix A and the matrix B are iteratively updated for several rounds, where any round of iterative update includes:

[0011] The first stream takes the upper left corner sub-matrix B of the unupdated part of matrix B 11 Send to the host;

[0012] The second stream will be the upper left corner sub-matrix A of the unupdated part of matrix A 11 Send to the host side, in parallel, the host side sends 11 Perform Cholesky decomposition to obtain the lower triangular matrix L B ;

[0013] The first stream receives L from the host B And store it in B 11 , in parallel, the host side according to L B To A 11 Normalize and get the matrix C A ;

[0014] The second stream receives C from the host side A And store it in A 11 , in parallel, the host side broadcasts L B ;

[0015] The first stream computes the matrix B 21 And send it to the host side matrix W B21 , in parallel, the second stream is based on Calculate the matrix A 21 , in parallel, the host side broadcasts C A , where W B21 For B 21 Auxiliary matrices of the same size, For L B The inverse conjugate transposed matrix, A 21 A 11 The bar matrix on the lower side, B 21 For B 11 Bar matrix on the lower side;

[0016] The second stream receives C from the host side A , the second stream according to B 21and C A Update A 21 , the second stream will be A 21 Matrix W sent to the host A21 , the second stream updates A 21 And send it to the host side matrix W A21 , in parallel, the host side broadcasts W B21 , where W A21 For A 21 Auxiliary matrices of the same size;

[0017] The host will W B21 Send to the first stream, in parallel, the host side calculates W B21 The transposed matrix U B12 , the host side broadcasts U B12 ;

[0018] The first stream updates the matrix B 22 、W AB and U BA , in parallel, the host side broadcasts W A21 , wherein said B 22 For B 21 The square matrix on the right, the W AB and U BA is the auxiliary matrix;

[0019] The host will W A21 Send to the second stream, in parallel, the host side calculates W A21 The transposed matrix U A12 , the host side broadcasts U A12 ;

[0020] The second stream updates the matrix W AB 、U BA 、A 22 and A 21 , wherein the A 22 A 21 The square on the right;

[0021] A 11 The rows and columns contained in A are marked as updated, and B 11 The included rows and columns are marked as updated in B.

[0022] In a possible implementation manner, the host side is based on L B To A 11 Normalize and get the matrix C A , including: according to the formula Calculate C A ,in, For L B The inverse matrix of .

[0023] In a possible implementation, the first stream calculates the matrix B 21 And send it to the host side matrix W B21 ,include:

[0024] The first stream receives L from the host B , the first stream according to the matrix Calculate matrix B 21 , the first-rate will B 21 Matrix W sent to the host B21 .

[0025] In a possible implementation, the second stream updates A 21 And send it to the host side matrix W A21 ,include:

[0026] The second stream receives C from the host side A , the second stream according to B 21 and C A Update A 21 , the second stream will be A 21 Matrix W sent to the host A21 .

[0027] In a possible implementation manner, the first stream updates the matrix B 22 、W AB and U BA ,include:

[0028] The first stream receives U B12 , the first stream according to W B21 and U B12 Update Matrix 1 22 , the first-rate will W B21 Fill in the matrix W AB The corresponding position, the first stream will U B12 Fill in the matrix U BA The corresponding position, where W AB With W B21 The number of rows is the same, the number of columns is W B21 2 times, U BA with U B12 The number of columns is the same, and the number of rows is U B12 2 times;

[0029] The second stream updates the matrix W AB 、U BA 、A 22 and A 21 ,include:

[0030] The second stream receives U A12 , the second stream will W A21Fill in the matrix W AB The corresponding position, the second flow will U A12 Fill in the matrix U BA The corresponding position, the second flow according to W AB and U BA Update matrix A 22 , the second stream according to B 21 and C A Update A again 21 .

[0031] In a possible implementation, the first stream is based on the matrix Calculate matrix B 21 , including: according to the formula Calculate B 21 .

[0032] In one possible implementation, the second stream is based on Calculate the matrix A 21 , including: according to the formula Calculate A 21 .

[0033] In one possible embodiment, the second stream is based on B 21 and C A Update A 21 , including: according to the formula Update A 21 .

[0034] In one possible implementation, the first stream is based on W B21 and U B12 Update matrix B 22 , including: According to formula B 22 =W B21 U B12 Update B 22 The second stream is based on W AB and U BA Update matrix A 22 , including: According to formula A 22 =W AB U BA Update A 22 The second stream is based on B 21 and C A Update A again 21 , including: according to the formula Update A again 21 .

[0035] In a possible implementation, the host end is a two-dimensional array of central processing units (CPUs); and the device end is a cluster of graphics processing units (GPUs).

[0036] This paper proposes a method for normalizing the generalized Hermitian matrix eigenvalue problem, which implements four types of parallel processing: 1. Parallel processing of host-side computation and data transfer between the host and the device. 2. Parallel processing of host-side communication and data transfer between the host and the device. 3. Parallel processing of device-side computation and data transfer between the host and the device. 4. Parallel processing of host-side communication and device-side computation. This parallel processing improves execution efficiency by approximately 30% compared to traditional normalization techniques. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 A flowchart of a method for normalizing the generalized Hermitian matrix eigenvalue problem disclosed in an embodiment of the present invention;

[0038] Figure 2 A flowchart of one round of iteration of a method for normalizing the generalized Hermitian matrix eigenvalue problem disclosed in an embodiment of the present invention;

[0039] Figure 3 A schematic diagram of an iterative update matrix disclosed in an embodiment of the present invention;

[0040] Figure 4 This is a performance comparison chart of the embodiment of the present invention in scenarios with different numbers of graphics processors;

[0041] Figure 5 This is a performance comparison chart for scenarios with different matrix sizes disclosed in an embodiment of the present invention. DETAILED DESCRIPTION

[0042] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only 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 making creative efforts shall fall within the scope of protection of the present invention.

[0043] Before introducing the embodiments of the present application in detail, the terms and symbols involved in the embodiments of the present application are first explained.

[0044] A Hermitian matrix, also known as a self-adjoint matrix, is a conjugate symmetric square matrix. In a Hermitian matrix, every element in the i-th row and j-th column is the complex conjugate of the element in the j-th row and i-th column.

[0045] Cholesky decomposition: refers to the decomposition of a positive Hermitian matrix into the product of a lower triangular matrix and its conjugate transposed matrix.

[0046] The traditional parallel generalized Hermitian normalization technique uses the method of first calculating the Cholesky decomposition of the positive definite matrix B and then calculating the matrix However, for the communication in the parallel solution of the Cholesky decomposition of the matrix B, for the solution of the matrix It is also necessary, which leads to repeated communication about matrix B. This patent proposes a new parallel heterogeneous block method for the normalization of the generalized Hermitian matrix characteristic problem. This method combines the Cholesky decomposition of matrix B with the solution of matrix In heterogeneous computing, because communication between devices cannot be carried out directly, data must first be transferred from the device to the host, then transmitted back to the device through communication on the host. Therefore, the standardized method proposed in this patent can reduce communication between hosts and between hosts and devices.

[0047] In parallel heterogeneous computing, in addition to avoiding unnecessary communication between hosts and between hosts and devices, the standardized technology of combining these two processes can also bring many benefits. By using two streams on the device side to control the Cholesky decomposition of matrix B and the matrix This allows the device to do work on matrix B while the host is doing work on matrix B. Related calculations; do matrix on the host side While the host is performing work related to , the device can be performing work related to matrix B, thus achieving asynchronous operation between the host and device. Specifically, four forms of host-device overlap can be achieved: overlap between host-side computation and data transfer between the host and device; overlap between host-side communication and data transfer between the host and device; overlap between device-side computation and data transfer between the host and device; and overlap between host-side communication and device-side computation.

[0048] In addition, this parallel heterogeneous normalization method merges multiple small-scale data processing steps into a larger-scale data processing, thereby increasing the parallelism and task load balancing of the method.

[0049] To facilitate understanding of the embodiments of the present invention, specific embodiments will be further explained below with reference to the accompanying drawings. The embodiments do not limit the embodiments of the present invention.

[0050] Figure 1 This is a flow chart of a method for normalizing the generalized Hermitian matrix characteristic problem disclosed in an embodiment of the present invention. Figure 1 As shown, the method is applied to a computing system, the computing system including a host side and a device side, the device side including a first stream and a second stream. In some possible implementations, the host side is a two-dimensional array of central processing units (CPUs), and the device side is a cluster of graphics processing units (GPUs).

[0051] In step 101 , a first stream receives a matrix B, and a second stream receives a matrix A, wherein the matrix A is a Hermitian matrix, the matrix B is a positive definite Hermitian matrix, and the sizes of the matrix A and the matrix B are n.

[0052] In step 102, the matrix A and the matrix B are updated by nt rounds of iterations, wherein: is the rounding up of x, and nb is the size of the algorithm block.

[0053] In some possible implementations, the size of the algorithm block is the same as the size of the CPU two-dimensional array storage block.

[0054] The steps of any round of iterative update in step 102 are as follows: Figure 2 shown.

[0055] Figure 2 This is a flowchart of one round of iteration of a method for normalizing the generalized Hermitian matrix characteristic problem disclosed in an embodiment of the present invention. Figure 2 As shown, stream1 represents the first stream on the device, and stream2 represents the second stream on the device. Steps marked with (CC-*) represent communication operations between the host and the device, steps marked with (H2D-*) represent host-to-device data transfers, steps marked with (D2H-*) represent device-to-host data transfers, and steps marked with (D-*) represent device-to-device operations. Steps followed by an asterisk (*) represent substeps of the Cholesky decomposition of matrix B.

[0056] For a variable that will be transmitted between the host side and the device side, in order to distinguish the variables between the host side and the device side in terms of expression, the variable prefixed with "d" represents the variable on the device side, and the variable without the "d" prefix represents the variable on the host side. It should be noted that for the same variable, having a "d" prefix and not having a "d" prefix only represent that it is located at different ends, but the variables represented by the two are essentially the same. For example, for matrix B and matrix dB, the former represents matrix B located at the host side, and the latter represents matrix B located at the device side. Both represent matrix B, and the "d" prefix is only set to facilitate the distinction between different ends in terms of expression during the communication process. In some possible implementations, since the programs on the host side and the device side are stored on different storage devices, they can use variable B to represent matrix B on their respective storage devices without causing conflicts.

[0057] It should be noted that within any of the following steps, the steps before and after the phrase "in parallel" are processed in parallel, while the steps on the same side of the phrase "in parallel" are processed in serial. For example, in the flow: "Step 1, Step 2, in parallel, Step 3," "Step 1, Step 2," and Step 3 are processed in parallel, meaning that Step 1 and Step 3 can be executed simultaneously, and Step 2 and Step 3 can also be executed simultaneously; whereas, Step 1 and Step 2 are processed in serial, meaning that Step 2 can only be executed after Step 1 is completed.

[0058] like Figure 2 As shown:

[0059] In step D2H-1, the first stream converts the upper left corner sub-matrix B of the unupdated part of matrix B into 11 Sent to the host, where B 11 The size of nb is nb.

[0060] In step D2H-2, the second stream converts the upper left corner sub-matrix A of the unupdated part of matrix A into 11 Send to the host side, in parallel, in step (1)*, the host side B 11 Perform Cholesky decomposition to obtain the lower triangular matrix L B , where A 11 The size of nb is nb.

[0061] In some possible implementations, in step (1)*, the host side sends a 11 Perform Cholesky decomposition by calling the corresponding function in the general linear algebra function library.

[0062] In step H2D-3, the first stream receives the calculated L from the host. B And store in dB 11, in parallel, in step (2), the host side according to L B To A 11 Normalize and get the matrix C A .

[0063] Specifically, in step (2), according to the formula Calculate C A ,in, For L B The inverse matrix of For L B The inverse conjugate transposed matrix of .

[0064] In some possible implementations, the host side performs the following operations according to L B To A 11 Normalization is performed by calling the corresponding function in the general linear algebra function library.

[0065] In step H2D-4, the second stream receives the calculated C from the host. A And store it in dA 11 , in parallel, in step CC-1, the host side broadcasts L B .

[0066] In step H2D-5, the first stream receives L from the host side. B And store it in dL B , in step (3)*, the first stream is calculated according to the matrix Calculate the dB matrix 21 , in step D2H-6, the first stream will be dB 21 Matrix W sent to the host B21 , in parallel, in step (4), the second stream is based on Calculate the matrix dA 21 , in parallel, in step CC-2, the host side broadcasts C A ,in, For L B The inverse conjugate transposed matrix, A 21 A 11 The bar matrix on the lower side, B 21 For B 11 The lower bar matrix, W B21 For B 21 Auxiliary matrices of the same size.

[0067] Specifically, in step (3)*, the first stream is calculated according to the formula Calculating dB 21 ; In step (4), the second stream is calculated according to the formula Calculate dA 21 .

[0068] In step H2D-7, the second stream receives C from the host. A and stored in dC A , in step (5), the second stream is based on B 21 and C A Update dA 21 In step D2H-8, the second stream converts dA 21 Matrix W sent to the host A21 , in parallel, in step CC-3, the host side broadcasts W B21 , where W A21 For A 21 Auxiliary matrices of the same size.

[0069] Specifically, in step (5), the second stream is calculated according to the formula Update dA 21 .

[0070] In step H2D-9, the host will B21 Send to the first stream, the first stream will W B21 Store to dW B21 , in parallel, in step CC-4, the host side calculates W B21 The transposed matrix U B12 In step CC-5, the host sends U B12 .

[0071] In step H2D-10, the first stream receives U from the host. B12 And store it in dU B12 , in step (6)*, the first stream is calculated according to dW B21 and dU B12 Update Matrix dB 22 In step D-1, the first stream will be dW B21 Fill in the matrix dW AB At the corresponding position, in step D-2, the first stream will B12 Fill in the matrix dU BA At the corresponding position, in parallel, in step CC-6, the host side broadcasts W A21 , where dW AB and dU BA is the auxiliary matrix, dW AB With dW B21 The number of rows is the same, the number of columns is dW B21 2 times, dU BA With dU B12 The number of columns is the same, the number of rows is dU B12 2 times, dB 22 dB 21 The square on the right.

[0072] Specifically, in step (6)*, the first stream is calculated according to formula B. 22 =W B21 U B12 Update dB 22 ; In step D-1, the first stream will dW B21 Fill in the matrix dW AB The right half; in step D-2, the first stream will dU B12 Fill in the matrix dU BA The upper part of the .

[0073] In step H2D-11, the host will A21 Send to the second stream, the second stream will W A21 Store to dW A21 , in parallel, in step CC-7, the host side calculates W A21 The transposed matrix U A12 In step CC-8, the host sends U A12 .

[0074] In step H2D-12, the second stream receives U from the host side. A12 And store it in dU A12 , in step D-3, the second stream will dW A21 Fill in the matrix dW AB At the corresponding position, in step D-4, the second stream will A12 Fill in the matrix dU BA At the corresponding position, in step (7), the second stream is based on W AB and U BA Update matrix A 22 , in step (8), the second stream is based on B 21 and C A Update A again 21 , where A 22 A 21 The square on the right.

[0075] Specifically, in step D-3, the second stream will dW A21 Fill in the matrix dW AB The left half; in step D-4, the second stream will dU A12 Fill in the matrix dU BA The lower half; in step (7), the second stream according to formula A C2 =W AB U BA Update dA 22 ; In step (8), the second stream is calculated according to the formula Update dA again 21 .

[0076] Finally, A 11The rows and columns contained in A are marked as updated, and B 11 The contained rows and columns are marked as updated in B, and the next iteration is performed.

[0077] In each iteration, the matrix A 11 、A 21 、A 22 The relative position relationship with A, matrix B 11 、B 21 、B 22 The relative position relationship between B and Figure 3 shown.

[0078] In this embodiment, steps (1)*, (3)*, and (6)* are sub-steps of Cholesky decomposition of matrix B; steps (2), (4), (5), (7), and (8) are sub-steps of normalizing matrix A to sub-steps.

[0079] In another embodiment, several rounds of iterations of step 102 can be described by the following pseudo code, wherein the step number in front of each line of code represents the step number in the Figure 2 The corresponding steps in ; variables prefixed with "d" represent variables on the device side, and variables without the "d" prefix represent variables on the host side; (stream1) represents the code executed in the first stream, and (stream2) represents the code executed in the second stream; Cbroadcast() represents the broadcast in the column communication domain, Rbroadcast() represents the broadcast in the row communication domain, and TransM() represents the matrix transposition operation between the row and column communication domains; A←B represents the transfer of matrix B to matrix A, and A→B represents the transfer of matrix A to matrix B; is rounded up to the nearest integer, nb is the size of the algorithm block, and n is the size of matrices A and B; CS is the set of column numbers, and RS is the set of row numbers, both initially {0, n-1}; StreamSynchronize(stream) means waiting for stream execution to complete before executing subsequent code; * in steps (D-1) to (D-4) represents a placeholder, indicating that the corresponding part of the matrix on the left side of the equal sign will not be modified. For example, step (D-1) represents replacing dW B Fill in dW AB The right half of dW AB The left half is not modified.

[0080]

[0081]

[0082] As shown in the above pseudocode, the present invention proposes a parallel hybrid method for converting the generalized Hermitian eigenvalue problem into the standard eigenvalue problem by combining Cholesky decomposition. It should be noted that in the two-dimensional block cyclic distribution of the matrix, each row and each column block of the matrix is only stored in a part of the processor, that is, stored in a row processor and a column processor respectively. This makes the increase of row and column blocks in the parallel method much more complicated than that of the serial method. The above pseudocode realizes the change of rows and columns of the local matrix distributed in each processor by continuously modifying the row set RS and column set CS so that they correspond to the unprocessed part of the local matrix.

[0083] When executing the program corresponding to the above pseudo code, two streams can be enabled at the same time on the device side, responsible for the Cholesky decomposition of matrix B and the matrix Calculation of . Figure 2 As shown, stream1 on the left is responsible for the calculation related to the Cholesky decomposition of matrix B, and stream2 on the right is responsible for the calculation related to the matrix The host is responsible for the relevant communication operations. Therefore, the operations related to the Cholesky decomposition of matrix B on the device side can be performed with the host side and the matrix Related operations or matrices The data transfer between the host and the device is carried out in parallel. The related operations can be performed in parallel with the host-side operations on matrix B or the host-device data transmission of matrix B. However, in the traditional method, the Cholesky decomposition of matrix B is first calculated, and then the matrix B is calculated. These processes can only be performed in serial mode. The method of the present invention realizes asynchronous cooperation between the device side and the host side, and realizes four forms of overlap (parallel):

[0084] 1. Overlap between host-side computation and host-device data transfer

[0085] When computing the Cholesky decomposition of the bt-th diagonal block of B and the normalization of the diagonal block of A, the diagonal blocks of B and A need to be transferred from the device to the host. After completing these two calculations on the host, the data needs to be transferred back to the device. By using asynchronous data transfer operations, the host calculations can be performed in parallel with the data movement between the host and the device. Figure 2 As shown, in step (1)*, the Cholesky decomposition of the diagonal blocks of matrix B can overlap with the device-to-host data transfer of the diagonal blocks in step (D2H-2). Similarly, the normalization process of the diagonal blocks of A in step (2) can overlap with the host-to-device data transfer of the diagonal blocks of B in step (H2D-3).

[0086] 2. Overlap between host-side communication and host-device data transmission

[0087] In the parallel heterogeneous generalized Hermitian eigenvalue problem standardization method, the communication between devices needs to be carried out through the MPI (Message Passing Interface) on the host side, so the data needs to be transmitted from the device side to the host side first, and then transmitted back to the device side after the communication is completed. The method described in the present invention combines the Cholesky decomposition of the matrix B and the matrix During the computational process, the host-side communication operations of matrix B can overlap with the host-to-device data transfer of matrix A. For example, the communication step (CC-1) of matrix B and the host-to-device data transfer step (H2D-4) of matrix A can overlap. Similarly, the host-side communication operations of matrix A can overlap with the host-to-device data movement of matrix B. For example, the communication step (CC-4) of matrix A and the host-to-device data transfer step (H2D-9) of matrix B can overlap.

[0088] 3. Overlap between device-side computing and host-device data transfer

[0089] like Figure 2 As shown, in step (D2H-6), the block matrix W of the matrix B is B The data transmission from the device side to the host side can overlap with the calculation of the partial matrix correction of matrix A in step (4).

[0090] 4. Overlap between host-side communication and device-side computation

[0091] The host-side communication operations for matrix B can overlap with the device-side calculations for matrix A. For example, the communication step (CC-3) for matrix B and the device-side calculation step (5) for matrix A can overlap. Similarly, the host-side communication operations for matrix A can overlap with the device-side calculations for matrix B. For example, the communication step (CC-2) for matrix A can overlap with the device-side calculation step (3)* for matrix B.

[0092] In a specific embodiment, a matrix of size 40000*40000 is normalized using the traditional method and the method of the present invention respectively. The test is performed under different numbers of GPUs and the time is recorded. The results are summarized as shown in Table 1.

[0093] Table 1 Comparison of test results improvement when using different numbers of GPUs for a 40,000-order matrix

[0094] Number of GPUs Traditional method (seconds) The method of the present invention (seconds) Performance improvements 4 151.07 120.25 20.40％ 8 92.25 65.40 29.11％ 16 54.56 38.48 29.47％ 32 43.93 30.62 30.29％

[0095] The bar chart corresponding to Table 1 is as follows Figure 4 As shown, traditional represents the results of the traditional method, corresponding to the dark gray bars on the left side of each group of bar graphs in the figure, and our represents the results of the method of the present invention, corresponding to the light gray bars on the right side of each group of bar graphs in the figure.

[0096] From Table 1 and Figure 4 It can be seen from the results that the method of the present invention has a performance improvement of 20%-30% compared to the traditional method.

[0097] In another specific embodiment, 32 GPUs are used to normalize matrices of different sizes using the traditional method and the method of the present invention, respectively. The test is performed and the time is recorded. The results are summarized as shown in Table 2.

[0098] Table 2 Comparison of test results improvement when using 32 GPUs for matrices of different sizes

[0099] Matrix size Traditional method (seconds) The method of the present invention (seconds) Performance improvements 40000 43.93 30.62 30.30％ 60000 113.13 76.14 32.70％ 80000 222.25 151.64 31.77％ 100000 365.50 259.87 28.90％

[0100] The bar chart corresponding to Table 2 is as follows Figure 5 As shown, traditional represents the results of the traditional method, corresponding to the dark gray bars on the left side of each group of bar graphs in the figure, and our represents the results of the method of the present invention, corresponding to the light gray bars on the right side of each group of bar graphs in the figure.

[0101] From Table 2 and Figure 5 It can be seen from the figure that the method of the present invention has a performance improvement of about 30% compared with the traditional method.

[0102] It should be noted that, in this document, relational terms such as first and second, etc., are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the terms "comprises," "comprising," or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, method, article, or device comprising a series of elements includes not only those elements, but also other elements not explicitly listed, or elements inherent to such process, method, article, or device. In the absence of further limitations, an element defined by the phrase "comprising a ..." does not exclude the presence of other identical elements in the process, method, article, or device comprising the element.

[0103] Those skilled in the art will understand that all or part of the steps to implement the above embodiments may be accomplished by hardware, or may be accomplished by a program instructing the relevant hardware, and the program may be stored in a computer-readable storage medium, and the above-mentioned storage medium may be a read-only memory, a disk, or an optical disk, etc.

[0104] The specific implementation methods described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific implementation method of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A generalized Hermitian matrix eigenvalue problem normalization method, applied to a computing system, the computing system comprising a host side and a device side, the device side comprising a first stream and a second stream, characterized in that: The method comprises: The first stream receives a matrix B, and the second stream receives a matrix A, wherein the matrix A is a Hermitian matrix and the matrix B is a positive definite Hermitian matrix; Perform several rounds of iterative updates on the matrix A and the matrix B, where any round of iterative updates includes: The first stream takes the upper left corner sub-matrix B of the unupdated part of matrix B 11 Send to the host; The second stream will be the upper left corner sub-matrix A of the unupdated part of matrix A 11 Send to the host side, in parallel, the host side sends 11 Perform Cholesky decomposition to obtain the lower triangular matrix L B ; The first stream receives L from the host B And store it in B 11 , in parallel, the host side according to L B To A 11 Normalize and get the matrix C A ; The second stream receives C from the host side A And store it in A 11 , in parallel, the host side broadcasts L B ; The first stream computes the matrix B 21 And send it to the host side matrix W B21 , in parallel, the second stream is based on Calculate the matrix A 21 , in parallel, the host side broadcasts C A , where W B21 For B 21 Auxiliary matrices of the same size, For L B The inverse conjugate transposed matrix, A 21 A 11 The bar matrix on the lower side, B 21 For B 11 Bar matrix on the lower side; The second stream receives C from the host side A , the second stream according to B 21 and C A Update A 21 , the second stream will be A 21 Matrix W sent to the host A21 , the second stream updates A 21 And send it to the host side matrix W A21 , in parallel, the host side broadcasts W B21 , where W A21 For A 21 Auxiliary matrices of the same size; The host will W B21 Send to the first stream, in parallel, the host side calculates W B21 The transposed matrix U B12 , the host side broadcasts U B12 ; The first stream updates the matrix B 22 、W AB and U BA , in parallel, the host side broadcasts W A21 , wherein said B 22 For B 21 The square matrix on the right, the W AB and U BA is the auxiliary matrix; The host will W A21 Send to the second stream, in parallel, the host side calculates W A21 The transposed matrix U C12 , the host side broadcasts U A12 ; The second stream updates the matrix W AB 、U BS 、A 22 and A 21 , wherein the A 22 A 21 The square on the right; A 11 The rows and columns contained in A are marked as updated, and B 11 The included rows and columns are marked as updated in B.

2. The method according to claim 1, characterized in that The host end is based on L B To A 11 Normalize and get the matrix C A , including: according to the formula Calculate C A ,in, For L B The inverse matrix of .

3. The method according to claim 1, characterized in that The first stream calculates the matrix B 21 And send it to the host side matrix W B21 ,include: The first stream receives L from the host B , the first stream according to the matrix Calculate matrix B 21 , the first-rate will B 21 Matrix W sent to the host B21 .

4. The method according to claim 1, wherein The second stream updates A 21 And send it to the host side matrix W A21 ,include: The second stream receives C from the host side A , the second stream according to B 21 and C A Update A 21 , the second stream will be A 21 Matrix W sent to the host A21 .

5. The method according to claim 1, wherein The first stream updates the matrix B 22 、W AB and U BA ,include: The first stream receives U B12 , the first stream according to W B21 and U B12 Update matrix B 22 , the first-rate will W B21 Fill in the matrix W AB The corresponding position, the first stream will U B12 Fill in the matrix U BA The corresponding position, where W AB With W B21 The number of rows is the same, the number of columns is W B21 2 times, U BA with U B12 The number of columns is the same, and the number of rows is U B12 2 times; The second stream updates the matrix W AB 、U BA 、A 22 and A 21 ,include: The second stream receives U A12 , the second stream will W A21 Fill in the matrix W AB The corresponding position, the second flow will U A12 Fill in the matrix U BA The corresponding position, the second flow according to W AB and U BA Update matrix A 22 , the second stream according to B 21 and C A Update A again 21 .

6. The method according to claim 3, characterized in that The first stream is based on the matrix Calculate matrix B 21 , including: according to the formula Calculate B 21 .

7. The method according to claim 1, characterized in that The second stream is based on Calculate the matrix A 21 , including: according to the formula Calculate A 21 .

8. The method according to claim 4, characterized in that The second stream is based on B 21 and C A Update A 21 , including: according to the formula Update A 21 .

9. The method according to claim 5, characterized in that The first stream is based on W B21 and U B12 Update matrix B 22 , including: According to formula B 22 =W B21 U B12 Update B 22 The second stream is based on W AB and U BA Update matrix A 22 , including: According to formula A 22 =W AB U BA Update A 22 The second stream is based on B 21 and C A Update A again 21 , including: according to the formula Update A again 21 .

10. The method according to claim 1, characterized in that The host side is a two-dimensional array of central processing units (CPUs); the device side is a cluster of graphics processing units (GPUs).

Images