A fast direct solution method for complex electromagnetic problems

Abstract

The invention discloses a heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems, comprising: an MPI process on a CPU side reads a model file and model parameters after three-dimensional geometric modeling of a target; all MPI processes are evenly divided into two groups, all basis functions obtained are divided into M groups, and M×M unit block matrices are generated by using the grouped basis functions to merge and form a block matrix Z; each MPI process obtains basis function information contained in the unit block matrix required by itself according to an acquisition rule related to the MPI process group; all MPI processes complete the LU decomposition process of the block matrix Z in a heterogeneous parallel decomposition method jointly used by a CPU side and a GPU side based on the basis function information obtained by themselves; a current coefficient matrix is ​​solved by back-substitution of a matrix equation group according to the decomposition result of the obtained block matrix Z; and the three-dimensional geometric model of the target is calculated according to the current coefficient matrix to obtain near-field data and far-field data.

Description

Technical Field

[0001] The present invention belongs to the field of computational electromagnetics, and in particular relates to a heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems. Background Art

[0002] The method of moments (MoM) in computational electromagnetics is a traditional numerical method with high computational accuracy and has been widely used to analyze complex electromagnetic radiation and scattering problems. Although the method of moments has high theoretical accuracy, a complex dense matrix of dimension N will be generated during the solution process. For the method of moments using direct solution, the storage complexity and computational complexity will increase exponentially with the increase of matrix dimension, which causes computer resources to be consumed quickly and the problem cannot be solved, which seriously limits the solution scale of the method of moments. Fortunately, with the development of high-performance parallel computing technology and supercomputers, distributed parallel computing has made the method of moments no longer limited to the computing resources of a single computer. By calling multiple computing cores, the demand for computing resources of the method of moments has been greatly alleviated, and the solution scale has been increased by several orders of magnitude, which makes the method of moments capable of solving some practical engineering problems with large electrical and multi-scale complex structures. However, it is obviously unrealistic to use thousands of processors for normalized calculations at will, so there is still a need for an efficient electromagnetic algorithm that can meet the computational accuracy with less computing resources.

[0003] Seeking accurate and efficient numerical modeling methods has always been an important topic of great concern in the field of computational electromagnetics. Until the turn of the century, the emergence of fast algorithms greatly reduced the storage complexity and computational complexity, making it possible to accurately handle large-scale and complex electromagnetic problems. Among them, the most representative fast algorithm is the Multilevel Fast Multipole Algorithm (MLFMA), which has the characteristics of high computational accuracy and fast computational speed, and is very suitable for radar cross section (RCS) analysis, target identification, stealth technology research and other fields. However, the fast algorithm using iterative solvers will inevitably face convergence problems when solving electromagnetic problems, especially when dealing with strong resonance, multi-scale and strong coupling boundary value problems with coexistence of micro and macro structures. The condition number of the matrix generated by the integral equation is usually poor, which leads to slow convergence or even non-convergence of the iterative solver. Therefore, before using the iterative solver for calculation, it is often necessary to select a suitable preconditioner technology for processing, and efficient preconditioner technology has always been one of the difficult problems that need to be solved internationally. In addition, the iterative solver is less efficient for problems with multiple right-hand sides (RHS), such as single-station RCS calculation, because each right-hand side in the RHS needs to be recalculated.

[0004] For the convergence and single-station RCS calculation problems encountered by the above iterative solver, the efficient fast direct solution algorithm may be a better choice. The fast direct solution method inverts the complete matrix, avoiding the convergence problem. At the same time, for multiple right-hand side problems, only back substitution is needed to obtain all solutions, which greatly improves the efficiency of solving the single-station RCS problem.

[0005] In addition, the rapid development of heterogeneous supercomputers in recent years has provided a broad space for the innovation and development of electromagnetic algorithms. With the help of heterogeneous computing, the solution scale of electromagnetic numerical algorithms can be further expanded under the premise of ensuring calculation accuracy. However, the advantages of heterogeneous computing cannot be directly converted into electromagnetic simulation capabilities. In order to fully utilize the computing resources of heterogeneous computers, it is necessary to develop corresponding heterogeneous parallel fast direct solution methods to realize the scattering and radiation calculations of various targets. Summary of the invention

[0006] In order to solve the above problems existing in the prior art, the present invention provides a heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems. The technical problem to be solved by the present invention is achieved by the following technical solutions:

[0007] A heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems is applied to a computing system including at least one CPU and a GPU, wherein the CPU as a CPU end contains multiple MPI processes; the computing system virtualizes the GPU in advance using Context technology so that one MPI process corresponds to one virtualized GPU end; the method comprises:

[0008] All the MPI processes on the CPU end read the model file and model parameters after the three-dimensional geometric modeling of the target is performed; wherein the model file includes the coordinate data after the three-dimensional geometric model of the target is meshed using triangular facets; the model parameters include the material parameters, boundary conditions and excitation sources of the target; wherein the number of MPI processes on the CPU end is an even number; a pair of triangular facets connected after meshing obtains an RWG basis function;

[0009] The CPU side evenly divides all MPI processes into a first group of processes and a second group of processes, divides all obtained basis functions into M groups, and uses a verification process between the grouped basis functions to generate M×M unit block matrices to merge and form a block matrix Z as a complex dense matrix; each MPI process obtains the basis function information contained in the unit block matrix required by itself according to the acquisition rules related to the MPI process group; wherein M is a natural number greater than 0;

[0010] All the MPI processes on the CPU side complete the LU decomposition process of the block matrix Z based on the basis function information obtained by themselves, using a heterogeneous parallel decomposition method jointly implemented by the CPU side and the GPU side, and obtain a decomposition result of the block matrix Z;

[0011] The CPU side solves the current coefficient matrix by back-substitution of the matrix equations according to the decomposition result of the block matrix Z;

[0012] The CPU calculates the three-dimensional geometric model of the target according to the current coefficient matrix to obtain the required near-field data and far-field data.

[0013] The present invention belongs to the field of computational electromagnetics and involves solving electromagnetic scattering and radiation problems of complex targets. In the solution provided by the embodiment of the present invention, a novel heterogeneous parallel fast direct solution algorithm based on ACA is proposed, and the CONTEXT technology provided by the GPU is used to propose a CPU / GPU heterogeneous parallel programming model based on MPI, which ensures the communication, calculation and storage balance between processes as much as possible;

[0014] The present invention has the following beneficial effects:

[0015] 1. Compared with the existing parallel moment method based on traditional LU decomposition, the present invention provides a heterogeneous parallel fast direct solution algorithm, which reduces the computational complexity and storage complexity required for parallel moment method calculation.

[0016] 2. Compared with the existing parallel fast method based on iterative solution, the heterogeneous parallel fast direct solution algorithm provided by the present invention can reduce the storage complexity and computational complexity of MoM, and avoid the problem of slow convergence or even non-convergence of the iterative solver.

[0017] 3. The present invention completes the filling, decomposition and compression of the matrix in an alternating manner, which not only supports the analysis of scattering problems of pure metals or pure media, but also supports the simulation analysis of electromagnetic scattering and radiation problems of complex targets (complex structures and complex media). The use of heterogeneous parallel fast direct solution methods can accelerate the solution process of electromagnetic scattering and radiation problems while using fewer computing resources. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 It is a schematic diagram of the load imbalance solution;

[0019] Figure 2 It is a schematic diagram of the load balancing solution;

[0020] Figure 3 A schematic flow chart of a heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems provided by an embodiment of the present invention;

[0021] Figure 4 A four-stage algorithm flow chart corresponding to a heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems provided by an embodiment of the present invention;

[0022] Figure 5 is a schematic diagram of a block matrix Z according to an embodiment of the present invention;

[0023] Figure 6a for Figure 5 The corresponding block matrix Z and process relationship diagram;

[0024] Figure 6b for Figure 6a Schematic diagram of the cell block matrix that needs to be filled in the 0 process;

[0025] Figure 6c for Figure 6a Schematic diagram of the unit block matrix that needs to be filled in process 2;

[0026] Figure 7 : is a mapping relationship between M×M unit block matrices and 6 processes in an embodiment of the present invention;

[0027] Figure 8A schematic diagram for understanding inter-process communication during the unit block matrix decomposition process of an embodiment of the present invention;

[0028] Fig. 9 Schematic diagram of a non-blocking sending process of process 1 and a non-blocking receiving process of process 2 in an embodiment of the present invention;

[0029] Fig.10 Schematic diagram of the second row 2 process execution unit block matrix decomposition processing process in an embodiment of the present invention;

[0030] Fig.11 A schematic diagram for understanding the subdivision form of matrix multiplication above the diagonal in an embodiment of the present invention;

[0031] Fig.12 It is a schematic diagram for understanding the task division of three-step GPU heterogeneous multiplication corresponding to the matrix multiplication subdivision form above the diagonal in an embodiment of the present invention;

[0032] Fig.13 Schematic diagram of the third row 2 process execution unit block matrix decomposition processing process in an embodiment of the present invention;

[0033] Fig.14 A schematic diagram for understanding the row-wise reorganization of a matrix in an embodiment of the present invention;

[0034] Fig.15 A schematic diagram for understanding the subdivision form of matrix multiplication below the diagonal line in an embodiment of the present invention;

[0035] Fig.16 It is a schematic diagram for understanding the task division of three-step GPU heterogeneous multiplication corresponding to the matrix multiplication subdivision form below the diagonal line in an embodiment of the present invention;

[0036] Fig.17 This is an example of the decomposition result of the block matrix Z in an embodiment of the present invention. DETAILED DESCRIPTION

[0037] The present invention is further described in detail below with reference to specific embodiments, but the embodiments of the present invention are not limited thereto.

[0038] In order to reduce the storage complexity and computational complexity of MoM, avoid the convergence problem of iterative solution algorithm, speed up the electromagnetic simulation speed of fine structure model and shorten the development cycle of civil equipment and major equipment, the present invention provides a heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems in view of the shortcomings of the prior art.

[0039] Before specifically describing the heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems, the binding relationship between the GPU (graphics processing unit) card and the process in the method is first described.

[0040] A heterogeneous node of a modern heterogeneous computer usually consists of multiple CPU (Central Processing Unit) cores and 1 to 4 GPU accelerator cards, so the number of GPU accelerator cards is often less than the number of CPUs. If an MPI process (MPI stands for Message Passing Interface) is bound to a GPU accelerator card, some processes will be unable to bind to the GPU accelerator card, which will cause load imbalance between processes.

[0041] In order to ensure load balancing between processes, the best way is to bind each process to the GPU accelerator card belonging to the process. In this case, Context technology is needed. This technology is mainly used to ensure load balancing of processes. It is usually used when the number of processes does not match the number of GPU cards. Figure 1 and Figure 2 Load imbalance and load balancing solutions are given respectively. Figure 1 It can be seen that some CPUs cannot match GPU cards. Figure 2 The second solution shown, Figure 2 As can be seen from the figure, by virtualizing the GPU card through the Context technology, each process can be matched to a virtual GPU card. The present invention regards each virtual GPU card as a GPU end.

[0042] Specifically, Context is similar to MPI process, managing the life cycle of resources allocated by the Driver layer. In other words, Context in GPU is equivalent to Program in CPU. There can be multiple Contexts on a GPU, and the Context addresses are isolated from each other. Therefore, Context technology can be used to virtualize a GPU card into multiple GPU cards, and then an MPI process is matched with a Context, so that each process can be guaranteed to have the same computing and storage capabilities as much as possible. The binding relationship between GPU card and process is as follows: Figure 2 shown.

[0043] Specifically, the heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems is applied to a computing system including at least one CPU and a GPU (i.e., a GPU accelerator card), wherein the CPU as a CPU end contains multiple MPI processes; the computing system virtualizes the GPU in advance using the Context technology, so that one MPI process corresponds to one virtualized GPU end;

[0044] The programming framework provided by the embodiment of the present invention is scalable and can be applied to multiple CPUs. For ease of understanding, the following description is based on a computing system including a group of CPUs and GPUs.

[0045] like Figure 3 As shown in FIG. 1 , the heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems may include the following steps S1 to S5. The implementation process can be divided into four parts, namely, reading the model (corresponding to S1), matrix filling and decomposition and compression (corresponding to S2 to S3), matrix back substitution (corresponding to S4) and post-processing calculation (corresponding to S5). The algorithm flow chart corresponding to the four stages is shown in FIG. Figure 4 As shown. Since the matrix filling, decomposition and compression in the second stage involve a large number of matrix multiplication operations, and this operation is the most time-consuming part in the program calculation process, it is very suitable to use heterogeneous computing. Therefore, the present invention mainly uses GPU components to accelerate the matrix multiplication operations in the second stage.

[0046] See also Figure 3 The specific steps of the method include:

[0047] S1, all MPI processes on the CPU end read the model file and model parameters after three-dimensional geometric modeling of the target;

[0048] The model file includes coordinate data of a three-dimensional geometric model of a target after meshing using triangular facets; the model parameters include material parameters, boundary conditions and excitation sources of the target; the number of MPI processes on the CPU side is an even number; a pair of triangular facets connected after meshing obtain an RWG basis function; the RWG basis function (RWG stands for Rao-Wilton-Glisson) is a basis function defined on adjacent planar triangular patches, and is a basis function commonly used in moment method processing of three-dimensional scattering problems.

[0049] The target of the embodiment of the present invention may include an object with a complex structure such as an antenna, and the material may be not only pure metal, but also pure dielectric, or a metal-medium mixed material. The present invention may add a plane wave excitation source, a coaxial wave port, or a rectangular wave port excitation source, and may calculate not only the scattering problem of a metal target, but also the radiation problem of a metal-medium mixed target.

[0050] The embodiment of the present invention can use LASPCEM simulation software to perform three-dimensional geometric modeling on the target. For the process of performing three-dimensional geometric modeling on the target and obtaining the model file and model parameters, please refer to the relevant technical understanding, which will not be described in detail here.

[0051] In the embodiment of the present invention, the model file and model parameters after the three-dimensional geometric modeling of the target are read by all MPI processes on the CPU side, and each process reads the same data.

[0052] S2, the CPU side evenly divides all MPI processes into a first group of processes and a second group of processes, divides all obtained basis functions into M groups, and generates M×M unit block matrices by using the inspection process between the grouped basis functions to merge and form a block matrix Z as a complex dense matrix; each MPI process obtains the basis function information contained in the unit block matrix required by itself according to the acquisition rules related to the MPI process group; wherein M is a natural number greater than 0;

[0053] Since the number of MPI processes on the CPU side is an even number, these MPI processes can be evenly divided into two groups, namely the first group of processes and the second group of processes, providing a new parallel strategy.

[0054] For ease of understanding, the following description is made by taking the CPU end containing 6 MPI processes as an example, but it should be known that the number of MPI processes on a CPU end in the embodiment of the present invention is not limited thereto.

[0055] The six MPI processes on the CPU side can be respectively recorded as process 0, process 1, process 2, process 3, process 4, and process 5. Process 0, process 1, and process 2 are classified into the first group of processes, and process 3, process 4, and process 5 are classified into the second group of processes.

[0056] Through S1, the CPU can obtain the RWG basis functions (hereinafter referred to as basis functions) obtained from all triangular facets, and then divide all the obtained basis functions into M groups according to the preset number of basis functions in each group, wherein the preset number of basis functions in each group can be any value between 500 and 1000, for example, if it is 500, then all basis functions are divided into M groups according to one group containing 500 basis functions.

[0057] The CPU side generates M×M unit block matrices by using the self-action and interaction test process between the grouped basis functions. The M×M unit block matrices can be combined to form a block matrix Z, which is a complex dense matrix in the moment method.

[0058] See also Figure 5After understanding the block matrix Z and the grouping of basis functions, the complete Z matrix can be pre-divided into multiple small matrices, namely unit block matrices. Assuming that the unit block matrix with row 3 and column 4 needs to be filled, it can be obtained by simply checking the basis functions of the third group with the basis functions of the fourth group.

[0059] For the verification process of the self-action and interaction between basis functions, please refer to the relevant technology and will not be explained here.

[0060] Each MPI process obtains the basis function information contained in the unit block matrix required by itself according to the acquisition rules related to the MPI process group, in preparation for filling the corresponding unit block matrix.

[0061] Each MPI process obtains the basis function information contained in the unit block matrix required by itself according to the acquisition rules related to the MPI process group, including:

[0062] In the block matrix Z, the basis function information of the unit block matrices on the diagonal and above the diagonal is mapped to each process in the first group of processes in sequence by column, and the basis function information of the unit block matrices on the diagonal and below the diagonal is mapped to each process in the second group of processes in sequence by row; wherein the diagonal refers to the line connecting the upper left corner and the lower right corner in the block matrix Z.

[0063] For ease of understanding, please refer to Figure 6a to Figure 6c Here, four processes, process 0, process 1, process 2, and process 3 are taken as examples. Process 0 and process 1 belong to the first group of processes, and process 2 and process 3 belong to the second group of processes. The above mapping process is explained by taking M×M unit block matrices as a 6*6 unit block matrix as an example.

[0064] Figure 6a yes Figure 5 The corresponding block matrix Z, the numbers on the top and left indicate the corresponding processes, Figure 6b Indicates the unit block matrix that needs to be filled in the 0 process, that is, the basis function information of which unit block matrices are mapped; Figure 6c It indicates the unit block matrices that need to be filled in process 2, and the unit block matrices that need to be filled in the remaining processes are no longer shown one by one.

[0065] See next Figure 7 Understand the mapping relationship between the M×M unit block matrices and the 6 processes, and each process obtains the basis function information of the corresponding unit block matrix. Figure 7 Here, 0, 1, 2, 3, 4, and 5 represent process 0, process 1, process 2, process 3, process 4, and process 5 respectively.

[0066] The present invention takes into account that data-dependent processes need to minimize the amount of data sent and the number of times the data is sent. In order to reduce communication between processes, the unit block matrices on the diagonal are mapped twice. Specifically, the basis function information of the two unit block matrices on the upper right and lower left of the diagonal is the same. If the calculation of the decomposition processing is performed by mapping once, the communication volume will be very large. Therefore, the unit block matrices on the diagonal are mapped twice, and one calculation is converted into two calculations. This can reduce communication, so that the same group of processes can use the first type of communication within the group to complete the sending of this part of information without sending it to another group of processes.

[0067] For ease of presentation, Figure 7 In the diagram, the same processes are represented by the same color.

[0068] S3, all the MPI processes on the CPU end use a heterogeneous parallel decomposition method jointly implemented by the CPU end and the GPU end based on the basis function information obtained by themselves to complete the LU decomposition process of the block matrix Z and obtain the decomposition result of the block matrix Z;

[0069] S3 can include:

[0070] Each process in the first group of processes performs decomposition processing on the portion of the block matrix Z located at the upper right of the diagonal line row by row based on the basis function information obtained by itself; wherein, in the upper right part, the decomposition processing of the unit block matrix at the diagonal position includes a filling link and a decomposition link, the decomposition processing of the unit block matrix at the non-diagonal position also includes a compression link after the decomposition link, and the decomposition processing of the unit block matrix other than the first row also includes an update link between the filling link and the decomposition link;

[0071] Each process in the second group of processes performs decomposition processing on the part of the block matrix Z located at the lower left of the diagonal line column by column based on the basis function information obtained by itself; wherein, in the lower left part, the decomposition processing of the unit block matrix at the diagonal position includes a filling link and a decomposition link, the decomposition processing of the unit block matrix at the non-diagonal position also includes a compression link after the decomposition link, and the decomposition processing of the unit block matrix other than the first column also includes an update link between the filling link and the decomposition link;

[0072] After the first group of processes completes the decomposition processing of the portion of the block matrix Z located at the upper right of the diagonal line row by row, and the second group of processes completes the decomposition processing of the portion of the block matrix Z located at the lower left of the diagonal line column by column, a decomposition result of the block matrix Z is obtained according to the decomposition processing results of all unit block matrices;

[0073] Among them, in the heterogeneous parallel decomposition method, the update phase of the unit block matrix is ​​executed by the GPU side, and the filling phase, decomposition phase and compression phase are executed by the CPU side. The above processes are all implemented by parallel programs

[0074] Combination Figure 8 Understand the above process and the inter-process communication in this process.

[0075] For ease of understanding, the decomposition process of the unit block matrix of the first row and the first column is first described.

[0076] (1) Decomposition of the unit block matrix of the first row includes:

[0077] In the first group of processes, the process corresponding to the unit block matrix at the diagonal position fills in the unit block matrix at the diagonal position using the basis function information obtained by itself, performs LU decomposition on the filled unit block matrix at the diagonal position, and obtains the corresponding L matrix and U matrix; the decomposed L matrix is ​​sent to the remaining processes in the first group of processes to complete the first type of communication of the first group of processes in the first row;

[0078] According to the order of the unit block matrices in the first row, the process corresponding to the unit block matrix at each non-diagonal position fills in the corresponding unit block matrix at the non-diagonal position using the basis function information obtained by itself, and decomposes the unit block matrix at the non-diagonal position filled by itself using the L matrix obtained by the first type of communication of the first group of processes in the first row to obtain the corresponding U matrix, and compresses the obtained U matrix; and if the unit block matrix at the non-diagonal position in the first row is located directly above a unit block matrix at a diagonal position in the block matrix Z, after compressing the obtained U matrix, the compression result is also sent to the matching process in the second group of processes to complete the second type of communication of the first group of processes in the first row;

[0079] For details, see Figure 8 In the first row of the upper right part, process 0 is the process corresponding to the unit block matrix at the diagonal position. Process 0 uses the basis function information obtained by itself to fill in the unit block matrix at the diagonal position, that is, the first unit block matrix in the first row. Next, the filled diagonal unit block matrix (that is, the first unit block matrix in the first row) is decomposed by LU to obtain the corresponding L matrix and U matrix; the calculation formula of LU decomposition is Z 11 =L 11 U 11 , Z 11 is the first row and first unit block matrix, Z 11 The subscripts of L represent the row and column numbers in the block matrix Z. The other matrices with subscripts in the embodiment of the present invention are similar and will not be described in detail below. 11 and U 11 It's Z 11The L matrix and U matrix are obtained by LU decomposition. For the specific process of LU decomposition, please refer to the relevant technical understanding and will not be described in detail here. Then the 0 process will decompose the L matrix 11 The matrix is ​​sent to processes 1 and 2 in the first group of processes to complete the first type of communication of the first group of processes in the first row. Figure 8 The black straight arrows from process 0 to process 1 and from process 0 to process 2 are indicated.

[0080] For process 1, the second unit block matrix of the first row is processed first. This unit block matrix is ​​the unit block matrix in the non-diagonal position. Process 1 uses the basis function information obtained by itself to fill in the second unit block matrix of the first row, and uses the L sent by the first type of communication of process 0 in the first row. 11 Matrix, decompose the first row and second unit block matrix filled in by process 1, and the calculation formula is U 12 =L 11 -1 Z 12 , where Z 12 Represents the first row and second unit block matrix, L 11 -1 Indicates L 11 The inverse matrix, U 12 Represents Z 12 The decomposed U matrix; since the second unit block matrix in the first row is a unit block matrix in a non-diagonal position, it is also necessary to get the U 12 The compression of the matrix is ​​performed by adopting the ACA (Adaptive Cross Approximation) algorithm in the embodiments of the present invention. Specifically, U 12 Compressed into the form of two U matrices multiplied to achieve dimensional change, U 12 Compressed into corresponding UuUv; For the second unit block matrix in the first row, since it is located directly above a unit block matrix at a diagonal position in the block matrix Z, combined with Figure 8 It can be seen that the second unit block matrix in the second row of the upper right part of the block matrix Z is a unit block matrix in a diagonal position. The second unit block matrix in the current first row and the unit block matrix in this diagonal position both belong to process 1 and are located directly above the unit block matrix in this diagonal position. Therefore, after compressing the U matrix obtained by the second unit block matrix in the first row, process 1 also sends the compression result to the second group of processes, that is, to processes 3, 4, and 5, to complete the second type of communication of the first group of processes in the first row; the information of the second type of communication is used by the receiving process to complete the update of the subsequent unit block matrix. For the second type of communication at this time, please refer to Figure 8 The green curved arrow from the second unit block matrix in the first row to the lower left area indicates.

[0081] For process 2, the third unit block matrix of the first row is processed first. This unit block matrix is ​​the unit block matrix of the non-diagonal position. Process 2 uses the basis function information obtained by itself to fill in the third unit block matrix of the first row, and uses the L sent by the first type of communication of process 0 in the first row. 11 Matrix, decompose the first row and third unit block matrix filled by process 2, and the calculation formula is U 13 =L 11 -1 Z 13 , where Z 13 Represents the third unit block matrix of the first row, U 13 Represents Z 13 The decomposed U matrix; since the third unit block matrix in the first row is a unit block matrix in a non-diagonal position, it is also necessary to get the U 13 The matrix is ​​compressed, specifically U 13 Compressed into the form of multiplication of two U matrices to achieve dimensional change, the compressed result is U 13 The corresponding UuUv, this compression result and U 12 The corresponding UuUv are different; for the third unit block matrix in the first row, since it is not directly above a unit block matrix in a diagonal position in the block matrix Z, after the two processes compress the U matrix obtained from the third unit block matrix in the first row, there is no need to send the compression result to the matching process in the second group of processes for the second type of communication.

[0082] Next, process 0 processes the fourth unit block matrix in the first row... and so on, until all the unit block matrices in the first row have completed the decomposition process.

[0083] It can be seen that in the first row, the decomposition process of the unit block matrix at the diagonal position includes a filling step and a decomposition step, and the first type of communication between the processes in the same group is performed after the decomposition, but the unit block matrix at the diagonal position does not have a compression step; the decomposition process of the unit block matrix at the non-diagonal position includes a filling step and a decomposition step, and a compression step is also included after the decomposition step, and if the unit block matrix at the non-diagonal position is located directly above a unit block matrix at a diagonal position in the block matrix Z, the compression result is also sent to the matching process in the second group of processes to complete the second type of communication between the first group of processes in the first row.

[0084] (2) Decomposition of the unit block matrix of the first column includes:

[0085] In the second group of processes, the process corresponding to the unit block matrix at the diagonal position fills in the unit block matrix at the diagonal position using the basis function information obtained by itself, performs LU decomposition on the filled unit block matrix at the diagonal position, and obtains the corresponding L matrix and U matrix; the decomposed U matrix is ​​sent to the remaining processes in the second group of processes to complete the first type of communication of the second group of processes in the first column;

[0086] According to the order of the unit block matrices in the first column, the process corresponding to the unit block matrix at each non-diagonal position fills in the corresponding unit block matrix at the non-diagonal position using the basis function information obtained by itself, and decomposes the unit block matrix at the non-diagonal position filled by itself using the U matrix obtained by the first type of communication of the second group of processes in the first column to obtain the corresponding L matrix, and compresses the obtained L matrix; and if the unit block matrix at the non-diagonal position in the first column is located to the left of a unit block matrix at a diagonal position in the block matrix Z, after compressing the obtained L matrix, the compression result is also sent to the matching process in the first group of processes to complete the second type of communication of the second group of processes in the first column.

[0087] For details, see Figure 8 In the first column of the lower left part, the 3 processes are the processes corresponding to the unit block matrix at the diagonal position. The 3 processes use the basis function information obtained by themselves to fill in the unit block matrix at the diagonal position, that is, the first unit block matrix in the first column. Next, the filled diagonal unit block matrix (that is, the first unit block matrix in the first column) is LU decomposed to obtain the corresponding L matrix and U matrix; where the calculation formula for LU decomposition is Z 11 =L 11 U 11 , Z 11 The first unit block matrix of the first column is the same as the first unit block matrix of the first row in the upper right part. 11 and U 11 It's Z 11 The L matrix and U matrix are obtained by LU decomposition. Then the 3rd process decomposes the U matrix 11 The matrix is ​​sent to processes 4 and 5 in the second group of processes to complete the first type of communication of the second group of processes in the first column. Figure 8 The black straight arrows from process 3 to process 4 and from process 3 to process 5 in the first column are used as illustrations.

[0088] For the 4 processes, the second unit block matrix in the first column is processed first. This unit block matrix is ​​the unit block matrix in the non-diagonal position. The 4 processes use the basis function information obtained by themselves to fill in the second unit block matrix in the first column, and use the U sent by the first type of communication of the 3 processes in the first column. 11Matrix, decompose the first column and second unit block matrix filled by 4 processes, the calculation formula is L 21 =Z 21 U 11 -1 , where Z 21 represents the second unit block matrix in the first column, U 11 -1 Indicates U 11 The inverse matrix, L 21 Represents Z 21 The decomposed L matrix; since the second unit block matrix in the first column is a unit block matrix in a non-diagonal position, it is also necessary to obtain the L 21 The matrix is ​​compressed, specifically L 21 Compressed into the form of multiplication of two L matrices to achieve dimensional change, L 21 Compressed into the corresponding LuLv; for the second unit block matrix in the first column, since it is located to the left of a unit block matrix at a diagonal position in the block matrix Z, combined with Figure 8 It can be seen that the second unit block matrix in the second column of the lower left part of the block matrix Z is a unit block matrix in a diagonal position. The second unit block matrix in the current first column and the unit block matrix in this diagonal position both belong to process 4 and are located directly to the left of the unit block matrix in this diagonal position. Therefore, after compressing the L matrix obtained by the second unit block matrix in the first column, process 4 also sends the compression result to the first group of processes, that is, to process 0, process 1, and process 2, to complete the second type of communication of the second group of processes in the first column; the information of the second type of communication is used by the receiving process to complete the subsequent update of the unit block matrix. For the second type of communication of the second group of processes in the first column, please refer to Figure 8 The red curved arrow in the first column indicates the direction from the second unit block matrix in the first column to the upper right area.

[0089] For the 5th process, the third unit block matrix in the first column is processed first. This unit block matrix is ​​the unit block matrix in the non-diagonal position. The 5th process uses the basis function information obtained by itself to fill in the third unit block matrix in the first column, and uses the U sent by the first type of communication of the 3rd process in the first column. 11 Matrix, decompose the first column and third unit block matrix filled by 5 processes, the calculation formula is L 31 =Z 31 U 11 -1 , where Z 31 Represents the third unit block matrix in the first column, L 31 Represents Z 31 The decomposed L matrix; since the third unit block matrix in the first column is a unit block matrix in a non-diagonal position, it is also necessary to obtain the L 31The matrix is ​​compressed, specifically L 31 Compressed into the form of multiplication of two L matrices to achieve dimensional change, the compressed result is L 31 The corresponding LuLv, this compression result and L 21 The corresponding LuLv is different; for the third unit block matrix in the first column, since it is not located to the left of a unit block matrix in a diagonal position in the block matrix Z, after the 5 processes compress the L matrix obtained by the third unit block matrix in the first column, there is no need to send the compression result to the matching process in the first group of processes by second type communication.

[0090] Next, three processes process the fourth unit block matrix in the first column... and so on, until all the unit block matrices in the first column have completed the decomposition process.

[0091] It can be seen that in the first column, the decomposition process of the unit block matrix at the diagonal position includes a filling step and a decomposition step, and the first type of communication between the processes in the same group is performed after the decomposition, but the unit block matrix at the diagonal position does not have a compression step; the decomposition process of the unit block matrix at the non-diagonal position includes a filling step and a decomposition step, and a compression step is also included after the decomposition step, and if the unit block matrix at the non-diagonal position is located to the left of a unit block matrix at a diagonal position in the block matrix Z, the compression result is also sent to the matching process in the first group of processes to complete the second type of communication between the second group of processes in the first column.

[0092] It should be noted that in the above decomposition process, after the program is started, each process can first fill in the first unit block matrix assigned to the process. After the unit block matrix at the diagonal position is decomposed, the process where the remaining unit block matrices in the first row are located uses the L matrix obtained by the first type of communication within the same process group to complete the row decomposition link, and completes the compression link for the decomposed L matrix. The process where the remaining unit block matrices in the first column are located uses the U matrix obtained by the first type of communication within the same process group to complete the column decomposition link, and completes the compression link for the decomposed U matrix. The first type of communication between processes is implemented through the broadcast of the message passing interface, because the remaining processes in the same process group can only perform row decomposition and column decomposition after obtaining the matrix sent by the process where the unit block matrix at the diagonal position is located. In order to overlap calculation and communication, the second type of communication uses non-blocking communication of the message passing interface.

[0093] The ACA algorithm used in the compression can approximate a low-rank matrix A of dimension n as the product of two small matrices of dimension n*k (k is much smaller than n), and the formula is A≈U*V. Both n and k are natural numbers greater than 0.

[0094] Next, the decomposition process of the unit block matrices of the remaining rows and columns except the unit block matrices of the first row and the first column will be described.

[0095] (I) Decomposition of the unit block matrix of each row except the first row, including:

[0096] The process corresponding to the unit block matrix at the diagonal position in the row fills in the unit block matrix at the diagonal position in the row using the basis function information obtained by itself, and updates the filled unit block matrix at the diagonal position in the row using the compression result sent by the second type of communication of the second group of processes and the compression result of all unit block matrices directly above the unit block matrix at the diagonal position in the row; performs LU decomposition on the updated unit block matrix at the diagonal position in the row to obtain the corresponding L matrix and U matrix; and sends the decomposed L matrix to the remaining processes in the first group of processes to complete the first type of communication of the first group of processes in the row;

[0097] According to the order of the unit block matrices of the row, the process corresponding to the unit block matrix at each non-diagonal position fills in the unit block matrices at all non-diagonal positions corresponding to the row using the basis function information obtained by itself, and updates the unit block matrices at all non-diagonal positions corresponding to the row filled in by the CPU end on the GPU end using the compression results of the unit block matrices at all non-diagonal positions corresponding to the row filled in and the unit block matrices directly above each, as well as the compression results sent by the second type of communication of the second group of processes before; decomposes the unit block matrices at the non-diagonal positions updated by itself using the L matrix obtained by the first type of communication of the first group of processes in the row to obtain the corresponding U matrix, and compresses the obtained U matrix; and, if the unit block matrix at the non-diagonal position of the row is located directly above a unit block matrix at a diagonal position in the block matrix Z, after compressing the obtained U matrix, the current compression result and the compression results of all the unit block matrices directly above are also sent to the matching process in the second group of processes to complete the second type of communication of the first group of processes in the row.

[0098] See also Figure 8 The following exemplary description is given in detail.

[0099] 1) For the second row

[0100] Process 1 uses the basis function information obtained by itself to fill in the unit block matrix in the diagonal position of the second row, that is, the second unit block matrix in the second row. It uses the compression result sent by the second type of communication of the second group of processes before, that is, the compression result LuLv sent by process 4, and the compression result UuUv of the second unit block matrix in the first row before process 1 to update the unit block matrix in the diagonal position of the filled second row, that is, the second unit block matrix in the second row. It is updated directly on the CPU side. The update formula is: Z after update 22 =Original Z 22 -[LuLv] 21 [UuUv] 12 Among them, the original Z 22 The second row and second unit block matrix is ​​filled out, [LuLv] 21 Indicates the compression result LuLv, [UuUv] sent by 4 processes through the second type of communication 12 It represents the compression result UuUv of the second unit block matrix of the first row before the 1st process. Then the updated second row second unit block matrix is ​​decomposed by LU, and decomposed into L 22 U 22 Since the second unit block matrix in the second row is a unit block matrix in the diagonal position, the first process needs to decompose the L 22 The matrix is ​​sent to process 0 and process 2, completing the first type of communication of the first group of processes in the second row. Figure 8 In order to overlap computation and communication, the second type of communication of process 1 in the first row and the first type of communication of process 1 in the second row are both implemented using non-blocking communication of the message passing interface. Therefore, processes other than process 0 and process 3 need to establish communication before filling the cell block matrix. The non-blocking sending process of process 1 and the non-blocking receiving process of process 2 are shown in Fig. 9 shown.

[0101] According to the order of the second row, the third unit block matrix of the second row needs to be processed. It is a unit block matrix in a non-diagonal position, corresponding to process 2. At this time, process 2 first uses the basis function information obtained by itself to fill in all the unit block matrices in the non-diagonal position corresponding to process 2 in the second row. The filling process is performed on the CPU side. Please refer to Fig.10 Example, 2 processes fill in all the cell blocks in the second row belonging to their own process matrix Z 23 , Z 26 , Z 29 , corresponding to the third unit block matrix in the second row, the sixth unit block matrix in the second row, and the ninth unit block matrix in the second row. After the filling is completed, the update phase is executed, which is performed on the GPU side. The update phase uses the filled Z 23 , Z26 , Z 29 Matrix, calculated Z 23 , Z 26 , Z 29 The compression result of the unit block matrix directly above the matrix is ​​[UuUv] 13 , [UuUv] 16 , [UuUv] 19 , it is also necessary to use the compression result sent by the second group of processes in the second type of communication, that is, the compression result sent by the 4 processes through the second type of communication [LuLv] 21 .

[0102] For the update process, for each row other than the first row, the process corresponding to each unit block matrix at each non-diagonal position fills in the unit block matrices at all non-diagonal positions corresponding to the row using the basis function information obtained by itself, and uses the compression results of the unit block matrices at all non-diagonal positions corresponding to the row filled in and the unit block matrices directly above each other, as well as the compression results sent by the second type of communication of the second group of processes, to update the unit block matrices at all non-diagonal positions corresponding to the row filled in by the CPU end on the GPU end, including:

[0103] For the process corresponding to the unit block matrix at the non-diagonal position in the row, the filled unit block matrices at all non-diagonal positions corresponding to the row are organized into Local_C data corresponding to the process by columns, the compression results of the unit block matrices directly above the unit block matrices at all non-diagonal positions corresponding to the row are organized into Local_B data corresponding to the process by columns, the compression results of the second type of communication sent by the second group of processes are organized into Local_A data corresponding to the process by columns, the Local_C data, Local_B data and Local_A data corresponding to the process are uploaded to the GPU end, the updated Local_C data is calculated on the GPU end based on the three-step GPU heterogeneous multiplication and downloaded to the CPU end, and the updated results of the unit block matrices at all non-diagonal positions corresponding to the process in the row are obtained after splitting by columns.

[0104] See also Fig.10 , taking 2 processes as an example to illustrate the updating process of the unit block matrix at the non-diagonal position, the corresponding Local_C data is composed of all the unit block matrices at the non-diagonal position filled by 2 processes in the second row, that is, the filled Z 23 , Z 26 , Z 29 The matrix is ​​organized into columns; Local_B data is Z 23 , Z 26 , Z 29The compression result of the unit block matrix directly above each matrix is ​​composed of columns, that is, [UuUv] 13 , [UuUv] 16 , [UuUv] 19 Composed by columns, for [UuUv] 13 , [UuUv] 16 , [UuUv] 19 Each of them will be split into a U u and a U v Then follow U u and U v The groups are arranged in columns; Local_A data is the compressed result sent by the second type of communication of the second group of processes [LuLv] 21 Composed of columns, Fig.10 Here [LuLv] 21 The two matrices in are represented as L u and L v . Fig.10 The bottom box shows the result of reorganizing the Local_A data, Local_B data, and Local_C data by columns. The purpose of reorganizing by columns is to give full play to the computing performance of the GPU and reduce the number of communications between the CPU and the GPU.

[0105] After the Local_C data, Local_B data and Local_A data corresponding to the 2nd process are uploaded to the GPU end, the GPU end calculates the updated Local_C data based on the three-step GPU heterogeneous multiplication. Fig.11 Schematic representation of the matrix multiplication subdivision above the diagonal, Fig.11 The part above the arrow is a schematic diagram of the matrix column reorganization process, and the part below the arrow is a three-step GPU heterogeneous multiplication process. In the three-step GPU heterogeneous multiplication, the first two steps are calculated using the cublaszgemm function that comes with the cuBLAS platform, and the third step uses a pre-written kernel function for heterogeneous calculation. By writing a kernel function and calling the cuBLAS math library provided by NVDIA, the matrix multiplication operation in the matrix decomposition process can be calculated, which can accelerate the MoM solution process and improve the computational performance of the fast and direct solution method.

[0106] For details, see Fig.11 , the three-step GPU heterogeneous multiplication includes the following steps:

[0107] Step 1: Use the cublaszgemm function in the cuBLAS platform to calculate Fig.11 Step ①L v ×U u Among them, L v From Local_A data, Uu Derived from Local_B data.

[0108] Step 2: Use the cublaszgemm function in the cuBLAS platform to calculate Fig.11 Step ② in L u ×L v U u .

[0109] Step 3, Fig.11 The updated Local_C in step ③ is Local_C = Local_C - Temp_C * U v It is calculated using kernel function. Fig.12 Represents the process of computing task allocation. Each matrix element in Local_C is mapped to a thread in the GPU. The computing process of each thread is i, j are the indices of each matrix element in the Local_C matrix. The Local_C_n(i,j) on the left side of the equal sign in the formula represents the updated result. Block 0 to Block n represent thread blocks in the GPU; Warp 0 to Warp m represent thread warps in the GPU.

[0110] The updated Local_C data is calculated through three-step GPU heterogeneous multiplication and downloaded to the CPU side. After being split by columns, the updated results of the unit block matrices of all non-diagonal positions corresponding to the 2 processes in the second row are obtained.

[0111] Then, the unit block matrices of all non-diagonal positions corresponding to the 2 processes in the second row are decomposed on the updated unit block matrices. The decomposition formula can be expressed as: U 2i =L 22 -1 Z 2i , where L 22 The matrix is ​​the L matrix decomposed by process 1 and sent to process 2 through the first type of communication of the first group of processes in the second row. Here Z 2i is the unit block matrix after the 2nd process update in the second row, U 2i is the corresponding decomposed U matrix, where i=3,6,9.... It can be understood that after the decomposition of all the unit block matrices in the non-diagonal positions corresponding to the 2 processes in the second row, the decomposed U matrix needs to be compressed, and for the third unit block matrix in the second row corresponding to the 2 processes, since it is located directly above a unit block matrix in a diagonal position, the compression result needs to be sent to the second group of processes using the second type of communication. At this time, what is sent via the second type of communication is the compression result of all the unit block matrices directly above the unit block matrix in the diagonal position of the third row, that is, the compression result of the third unit block matrix in the first row [UuUv] 13And the compression result of the second row and third unit block matrix [UuUv] 23 , they all belong to 2 processes.

[0112] The above is an explanation of the decomposition process of the second row using process 2 as an example. The remaining processes in the second row perform similar processing until all unit block matrices in the second row have completed the decomposition process.

[0113] 2) For the third row

[0114] See also Fig.13 , taking 2 processes as an example, the diagonal position unit block matrix is ​​updated on the CPU side. For details, please refer to the previous article for understanding. For the non-diagonal position unit block matrix, it is updated on the GPU side. At this time, 2 processes fill in the unit block matrix Z of the third row 36 , Z 39 , Z 3n , but the update needs to be performed twice, and the result of the first update is represented by Z 3i -[LuLv] 31 [UuUv] 1i , the second update result is represented by Z 3i -[LuLv] 32 [UuUv] 2i , where i = 6, 9, ... n; Z in the second update 3i Z after the first update 3i Each row update process is as follows Fig.11 and Fig.12 Similarly, after two update calculations are completed, the last updated unit block matrix Z is downloaded to the CPU side. For each updated Z, the decomposition process is first performed, and then the compression process is performed.

[0115] The decomposition and compression process is similar to the rest of the rows. In addition, the third unit block matrix in the third row, as the unit block matrix in the diagonal position, needs to complete the first type of communication with the first group of processes. The remaining unit block matrices in the third row, if they are located directly above the unit block matrix in the diagonal position, need to complete the second type of communication with the second group of processes. The specific description will not be repeated here.

[0116] (ii) Decomposition of the unit block matrix of each column except the first column includes:

[0117] The process corresponding to the unit block matrix at the diagonal position in the column fills in the unit block matrix at the diagonal position in the column using the basis function information obtained by itself, and updates the filled unit block matrix at the diagonal position in the column using the compression result sent by the second type of communication of the first group of processes and the compression result of all unit block matrices to the left of the unit block matrix at the diagonal position in the column; performs LU decomposition on the updated unit block matrix at the diagonal position in the column to obtain the corresponding L matrix and U matrix; and sends the decomposed U matrix to the remaining processes in the second group of processes to complete the first type of communication of the second group of processes in the column;

[0118] According to the order of the unit block matrices of the column, the process corresponding to the unit block matrix at each non-diagonal position fills in the unit block matrices at all non-diagonal positions corresponding to itself in the column using the basis function information obtained by itself, and updates the unit block matrices at all non-diagonal positions corresponding to itself in the row filled in by the CPU end on the GPU end using the compression results of the unit block matrices at all non-diagonal positions corresponding to itself in the filled column and the unit block matrices to the immediate left of each, as well as the compression results sent by the second type of communication of the first group of processes before; decomposes the unit block matrices at the non-diagonal positions updated by itself using the U matrix obtained by the first type of communication of the second group of processes in the column to obtain the corresponding L matrix, and compresses the obtained L matrix; and, if the unit block matrix at the non-diagonal position of the column is located to the immediate left of a unit block matrix at a diagonal position in the block matrix Z, after compressing the obtained L matrix, the current compression result and the compression results of all unit block matrices to the immediate left are also sent to the matching process in the first group of processes to complete the second type of communication of the second group of processes in the column.

[0119] See also Figure 8 The following exemplary description is given in detail.

[0120] 1) For the second column

[0121] Process 4 uses the basis function information obtained by itself to fill in the unit block matrix in the diagonal position of the second column, that is, the second unit block matrix in the second column. It uses the compression result sent by the second type of communication of the first group of processes, that is, the compression result UuUv sent by process 1, and the compression result LuLv of the second unit block matrix in the first column by process 4 before to update the unit block matrix in the diagonal position of the filled second column, that is, the second unit block matrix in the second column. Specifically, it is updated on the CPU side. The update formula is: Z after update 22 =Original Z 22 -[LuLv] 21 [UuUv] 12 Among them, the original Z 22The second row and second unit block matrix is ​​filled out, [UuUv] 12 Indicates the compression result UuUv, [LuLv] sent by process 1 through the second type of communication 21 It represents the compression result LuLv of the second unit block matrix of the first column before the 4th process. Then the updated second column second unit block matrix is ​​decomposed into L 22 U 22 Since the second unit block matrix in the second column is a unit block matrix in the diagonal position, the 4-process needs to decompose the U 22 The matrix is ​​sent to process 3 and process 5, completing the first type of communication of the second group of processes in the second column. At this time, process 4 completes the first type of communication of the second group of processes in the second column. Figure 8 As shown by the purple arrow in the second unit block matrix in the second column. In order to overlap the computation and communication, the second type of communication of the four processes in the first column and the first type of communication of the four processes in the second column are both implemented using non-blocking communication of the message passing interface.

[0122] According to the order of the second column, the third unit block matrix in the second column needs to be processed. It is a unit block matrix in a non-diagonal position, corresponding to 5 processes. At this time, 5 processes first use the basis function information obtained by themselves to fill in all the unit block matrices in the second column corresponding to 5 processes in non-diagonal positions. The filling link is performed on the CPU side. 5 processes fill in all the unit block matrices Z belonging to their own processes in the second column. 32 , Z 62 , Z 92 , corresponding to the third unit block matrix in the second column, the sixth unit block matrix in the second column, and the ninth unit block matrix in the second column. After the filling is completed, the update phase is executed, which is performed on the GPU side. The update phase uses the filled Z 32 , Z 62 , Z 92 Matrix, filled with Z 32 , Z 62 , Z 92 The compression result of the unit block matrix on the left side of the matrix is ​​[LuLv] 31 、[LuLv] 61 、[LuLv] 91 , it is also necessary to use the compression result sent by the second type of communication of the first group of processes, that is, the compression result sent by process 1 through the second type of communication [UuUv] 12 .

[0123] For the update process, the process corresponding to each unit block matrix at a non-diagonal position fills in the unit block matrices at all non-diagonal positions corresponding to itself in the column using the basis function information obtained by itself, and uses the compression results of the unit block matrices at all non-diagonal positions corresponding to itself in the filled column and the unit block matrices to the left of each other, as well as the compression results sent by the second type of communication of the first group of processes, to update the unit block matrices at all non-diagonal positions corresponding to itself in the row filled in by the CPU end on the GPU end, including:

[0124] For the process corresponding to the unit block matrix at the non-diagonal position in the column, the filled unit block matrices at all non-diagonal positions corresponding to the column are organized into rows to form Local_C data corresponding to the process, and the compression results of the unit block matrices to the left of all the unit block matrices at the non-diagonal positions corresponding to the column are organized into rows to form Local_B data corresponding to the process, and the compression results of the second type of communication sent by the second group of processes are organized into rows to form Local_A data corresponding to the process, and the Local_C data, Local_B data and Local_A data corresponding to the process are uploaded to the GPU end, and the updated Local_C data is calculated on the GPU end based on the three-step GPU heterogeneous multiplication and downloaded to the CPU end, and the updated results of the unit block matrices at all non-diagonal positions corresponding to the process in the column are obtained after splitting by rows.

[0125] For 5 processes, see Fig.14 The row-wise reorganization of the matrix shown in the figure on the left, and Fig.15 The matrix multiplication subdivision form of the matrix block below the diagonal is shown.

[0126] The corresponding Local_C data is composed of rows of the cell block matrix of all non-diagonal positions filled by the 5 processes in the second column, that is, the filled Z 32 , Z 62 , …, Z m2 The matrix is ​​organized in rows; Local_A data is Z 32 , Z 62 …, Z m2 The compression result of the unit block matrix to the left of each matrix is ​​composed of columns, that is, [LuLv] 31 、[LuLv] 61 , …, [LuLv] m1 By row, for [LuLv] 31 、[LuLv] 61 , …, [LuLv] m1 Each of them will be split into an L u and an L v Then follow L uand L v The groups are arranged in rows; Local_B data is the compressed result sent by the second type of communication of the first group of processes [UuUv] 12 The purpose of row-by-row reorganization is to give full play to the computing performance of GPU and reduce the number of communications between CPU and GPU.

[0127] After the Local_C data, Local_B data and Local_A data corresponding to the 5 process are uploaded to the GPU side, the GPU side calculates the updated Local_C data based on the three-step GPU heterogeneous multiplication. Fig.15 The matrix multiplication subdivision form below the diagonal line is shown. In the three-step GPU heterogeneous multiplication, the first two steps are calculated using the cublaszgemm function provided by the cuBLAS platform, and the third step uses a pre-written kernel function for heterogeneous calculation.

[0128] Specifically, this part of the three-step GPU heterogeneous multiplication includes the following steps:

[0129] Step 1: Use the cublaszgemm function in the cuBLAS platform to calculate Fig.15 Step ①L v ×U u Among them, L v From Local_A,U u From Local_B.

[0130] Step 2: Use the cublaszgemm function in the cuBLAS platform to calculate Fig.15 Step ② in L v U u ×U v .

[0131] Step 3, Fig.15 The updated Local_C in step ③ is Local_C = Local_C - L u *Temp_C is calculated by kernel function, and its tasks are divided as follows Fig.16 shown.

[0132] The updated Local_C data is calculated through three-step GPU heterogeneous multiplication and downloaded to the CPU side. After being split by rows, the updated results of the unit block matrices of all non-diagonal positions corresponding to the 5 processes in the second column are obtained.

[0133] Then, the unit block matrices of all non-diagonal positions corresponding to the 5 processes in the second column perform the decomposition step on the updated unit block matrix. The decomposition formula can be expressed as: L i2 =Z i2 U 22-1 , where U 22 The matrix is ​​the U matrix decomposed by process 4 and sent to process 5 through the first type of communication of the second group of processes in the second column. Here Z i2 is the unit block matrix after the 5-process update in the second column, L i2 is the corresponding decomposed U matrix, where i=3,6,9.... It can be understood that after the decomposition of all the unit block matrices in the non-diagonal positions corresponding to the 5 processes in the second column, the decomposed L matrix must be compressed, and for the third unit block matrix in the second column corresponding to the 5 processes, since it is located to the left of a unit block matrix in a diagonal position, the compression result needs to be sent to the first group of processes using the second type of communication. At this time, what is sent via the second type of communication is the compression result of all the unit block matrices to the left of the unit block matrix in the diagonal position of the third column, that is, the compression result of the third unit block matrix in the first column [LuLv] 31 And the compression result of the second column and the third unit block matrix [LuLv] 32 , they all belong to 5 processes.

[0134] The above is an explanation of the decomposition process of the second column using the 5th process as an example. The remaining processes in the second column are processed similarly until all the unit block matrices in the second column have completed the decomposition process.

[0135] 2) For the third column

[0136] See also Fig.14 In the right figure, taking 5 processes as an example, the diagonal unit block matrix is ​​updated on the CPU side. For details, please refer to the previous article. For the non-diagonal unit block matrix, it is updated on the GPU side. At this time, 5 processes fill in the non-diagonal unit block matrix Z in the third column. 63 , Z 93 , ..., Z m3 , the update also needs to be performed twice, and the result of the first update is represented by Z i3 -[LuLv] 13 [UuUv] i1 , the second update result is represented by Z i3 -[LuLv] 23 [UuUv] i2 , where i = 6, 9, ... n; Z in the second update i3 Z after the first update i3 Similarly, after the two update calculations are completed, the last updated unit block matrix Z is downloaded to the CPU. For each updated Z, the decomposition process is first performed, and then the compression process is performed.

[0137] The decomposition and compression process is similar to the rest of the rows, and the third unit block matrix in the third row, as the unit block matrix in the diagonal position, also needs to complete the first type of communication with the first group of processes. The remaining unit block matrices in the third row, if they are located directly above the unit block matrix in the diagonal position, also need to complete the second type of communication with the second group of processes.

[0138] Through the processing of S3, the decomposition result of the block matrix Z can be obtained, including the corresponding L matrix and U matrix.

[0139] It can be seen that the embodiment of the present invention realizes parallel heterogeneous solution in the update process of matrix multiplication. Specifically: the embodiment of the present invention provides two compressed matrix reorganization methods, one is the column-based reorganization form, and the other is the row-based reorganization form. The significance of matrix reorganization is to reduce the multiple communications between the CPU and the GPU to three communications, and at the same time synthesize the scattered matrices into a whole, which can give full play to the computing performance of the GPU. Moreover, the embodiment of the present invention provides a task division form of matrix multiplication, which subdivides the original update part into three steps of GPU heterogeneous multiplication after matrix reorganization. The first two steps are calculated using the cublaszgemm function that comes with the cuBLAS platform, and the third step uses the kernel function written by the present invention to perform heterogeneous calculations. The embodiment of the present invention also provides a kernel function execution strategy for GPU calculation matrix multiplication. As Fig.12 As shown, according to the matrix construction strategy, Temp_C*U v or L u *Temp_C is split into multiple groups of small matrix multiplications, and GPU heterogeneous computing is performed by mapping multiple small matrix multiplications to one or more thread blocks in the thread grid. The heterogeneous parallel fast direct solution algorithm provided by the present invention completes the filling, decomposition and compression of the matrix in an alternating manner, and can be used for electromagnetic scattering and radiation problem analysis of complex structures and complex media targets.

[0140] S4, the CPU side solves the current coefficient matrix by back-substitution of the matrix equations according to the decomposition result of the block matrix Z;

[0141] It can be understood by those skilled in the art that the essence of the moment method is to solve the linear equation ZI=V, where Z is a complex dense matrix of dimension N, I is the current coefficient matrix to be solved, and V is a known excitation matrix. The solution process of ZI=V can be further written as LUI=V, where LU is the L matrix and U matrix after the decomposition of the Z matrix. The present invention uses the ACA algorithm and heterogeneous computing to accelerate the solution process of linear equations.

[0142] S4 may include the following steps:

[0143] S41, the CPU obtains the L matrix and the U matrix decomposed from the block matrix Z from the decomposition result of the block matrix Z, and calculates the excitation matrix V according to the formula LUI=V based on the original solution formula ZI=V of the moment method; wherein I represents the current coefficient matrix to be solved;

[0144] in, Fig.17 is an example of the decomposition result of the block matrix Z. The upper right part is the decomposed U matrix, and the lower left part is the decomposed L matrix.

[0145] The calculation formula of the excitation matrix V includes:

[0146]

[0147] Among them, V J,m is the matrix element corresponding to the mth current basis function; V M,m is the matrix element corresponding to the mth magnetic flow basis function; f m is the RWG basis function; E inc is the incident electric field; H inc is the incident magnetic field; represents the equilateral triangle in the mth RWG basis function; Represents the negative triangle in the mth RWG basis function.

[0148] S42, setting UI as the intermediate process matrix Y, and calculating the intermediate process matrix Y according to the formula LY=V;

[0149] S43, calculating the current coefficient matrix I according to the formula UI=Y.

[0150] Through the above processing, the current coefficient matrix I for the target can be calculated.

[0151] S5, the CPU calculates the three-dimensional geometric model of the target according to the current coefficient matrix to obtain the required near-field data and far-field data.

[0152] According to the calculation results of S4, the target model is post-processed to calculate the required near-field and far-field data, and visual charts can be generated.

[0153] The present invention belongs to the field of computational electromagnetics and involves solving electromagnetic scattering and radiation problems of complex targets. In the solution provided by the embodiment of the present invention, a novel heterogeneous parallel fast direct solution algorithm based on ACA is proposed, and the CONTEXT technology provided by the GPU is used to propose a CPU / GPU heterogeneous parallel programming model based on MPI, which ensures the communication, calculation and storage balance between processes as much as possible;

[0154] The present invention has the following beneficial effects:

[0155] 1. Compared with the existing parallel moment method based on traditional LU decomposition, the present invention provides a heterogeneous parallel fast direct solution algorithm, which reduces the computational complexity and storage complexity required for parallel moment method calculation.

[0156] 2. Compared with the existing parallel fast method based on iterative solution, the heterogeneous parallel fast direct solution algorithm provided by the present invention can reduce the storage complexity and computational complexity of MoM, and avoid the problem of slow convergence or even non-convergence of the iterative solver.

[0157] 3. The present invention completes the filling, decomposition and compression of the matrix in an alternating manner, which not only supports the analysis of scattering problems of pure metals or pure media, but also supports the simulation analysis of electromagnetic scattering and radiation problems of complex targets (complex structures and complex media). The use of heterogeneous parallel fast direct solution methods can accelerate the solution process of electromagnetic scattering and radiation problems while using fewer computing resources.

[0158] The above description is only a preferred embodiment of the present invention and is not intended to limit the protection scope of the present invention. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention are included in the protection scope of the present invention.

Claims

1. A heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems, characterized in that: The method is applied to a computing system including at least one CPU and a GPU, wherein the CPU as a CPU end includes multiple MPI processes; the computing system virtualizes the GPU in advance using the Context technology so that one MPI process corresponds to one virtualized GPU end; the method includes: All the MPI processes on the CPU end read the model file and model parameters after the three-dimensional geometric modeling of the target is performed; wherein the model file includes the coordinate data after the three-dimensional geometric model of the target is meshed using triangular facets; the model parameters include the material parameters, boundary conditions and excitation sources of the target; wherein the number of MPI processes on the CPU end is an even number; a pair of triangular facets connected after meshing obtains an RWG basis function; The CPU side evenly divides all MPI processes into a first group of processes and a second group of processes, divides all obtained basis functions into M groups, and uses a verification process between the grouped basis functions to generate M×M unit block matrices to merge and form a block matrix Z as a complex dense matrix; each MPI process obtains the basis function information contained in the unit block matrix required by itself according to the acquisition rules related to the MPI process group; wherein M is a natural number greater than 0; All the MPI processes on the CPU side complete the LU decomposition process of the block matrix Z based on the basis function information obtained by themselves, using a heterogeneous parallel decomposition method jointly implemented by the CPU side and the GPU side, and obtain a decomposition result of the block matrix Z; The CPU side solves the current coefficient matrix by back-substitution of the matrix equations according to the decomposition result of the block matrix Z; The CPU calculates the three-dimensional geometric model of the target according to the current coefficient matrix to obtain the required near-field data and far-field data.

2. The heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems according to claim 1 is characterized in that: The acquisition rules related to the MPI process group include: In the block matrix Z, the basis function information of the unit block matrices on the diagonal and above the diagonal is mapped to each process in the first group of processes in sequence by column, and the basis function information of the unit block matrices on the diagonal and below the diagonal is mapped to each process in the second group of processes in sequence by row; wherein the diagonal refers to the line connecting the upper left corner and the lower right corner in the block matrix Z.

3. The heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems according to claim 1 is characterized in that: All the MPI processes on the CPU side complete the LU decomposition process of the block matrix Z based on the basis function information obtained by themselves, using a heterogeneous parallel decomposition method jointly implemented by the CPU side and the GPU side, and obtain the decomposition result of the block matrix Z, including: Each process in the first group of processes performs decomposition processing on the portion of the block matrix Z located at the upper right of the diagonal line row by row based on the basis function information obtained by itself; wherein, in the upper right part, the decomposition processing of the unit block matrix at the diagonal position includes a filling link and a decomposition link, the decomposition processing of the unit block matrix at the non-diagonal position also includes a compression link after the decomposition link, and the decomposition processing of the unit block matrix other than the first row also includes an update link between the filling link and the decomposition link; Each process in the second group of processes performs decomposition processing on the part of the block matrix Z located at the lower left of the diagonal line column by column based on the basis function information obtained by itself; wherein, in the lower left part, the decomposition processing of the unit block matrix at the diagonal position includes a filling link and a decomposition link, the decomposition processing of the unit block matrix at the non-diagonal position also includes a compression link after the decomposition link, and the decomposition processing of the unit block matrix other than the first column also includes an update link between the filling link and the decomposition link; After the first group of processes completes the decomposition processing of the portion of the block matrix Z located at the upper right of the diagonal line row by row, and the second group of processes completes the decomposition processing of the portion of the block matrix Z located at the lower left of the diagonal line column by column, a decomposition result of the block matrix Z is obtained according to the decomposition processing results of all unit block matrices; Among them, in the heterogeneous parallel decomposition method, the updating stage of the unit block matrix is ​​executed by the GPU side, and the filling stage, decomposition stage and compression stage are executed by the CPU side.

4. The heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems according to claim 3 is characterized in that: The decomposition process of the unit block matrix of the first row includes: In the first group of processes, the process corresponding to the unit block matrix at the diagonal position fills in the unit block matrix at the diagonal position using the basis function information obtained by itself, performs LU decomposition on the filled unit block matrix at the diagonal position, and obtains the corresponding L matrix and U matrix; the decomposed L matrix is ​​sent to the remaining processes in the first group of processes to complete the first type of communication of the first group of processes in the first row; According to the order of the unit block matrices in the first row, the process corresponding to the unit block matrix at each non-diagonal position fills in the corresponding unit block matrix at the non-diagonal position using the basis function information obtained by itself, and decomposes the unit block matrix at the non-diagonal position filled by itself using the L matrix obtained by the first type of communication of the first group of processes in the first row to obtain the corresponding U matrix, and compresses the obtained U matrix; and if the unit block matrix at the non-diagonal position in the first row is located directly above a unit block matrix at a diagonal position in the block matrix Z, after compressing the obtained U matrix, the compression result is also sent to the matching process in the second group of processes to complete the second type of communication of the first group of processes in the first row; The decomposition process of the unit block matrix of the first column includes: In the second group of processes, the process corresponding to the unit block matrix at the diagonal position fills in the unit block matrix at the diagonal position using the basis function information obtained by itself, performs LU decomposition on the filled unit block matrix at the diagonal position, and obtains the corresponding L matrix and U matrix; the decomposed U matrix is ​​sent to the remaining processes in the second group of processes to complete the first type of communication of the second group of processes in the first column; According to the order of the unit block matrices in the first column, the process corresponding to the unit block matrix at each non-diagonal position fills in the corresponding unit block matrix at the non-diagonal position using the basis function information obtained by itself, and decomposes the unit block matrix at the non-diagonal position filled by itself using the U matrix obtained by the first type of communication of the second group of processes in the first column to obtain the corresponding L matrix, and compresses the obtained L matrix; and if the unit block matrix at the non-diagonal position in the first column is located to the left of a unit block matrix at a diagonal position in the block matrix Z, after compressing the obtained L matrix, the compression result is also sent to the matching process in the first group of processes to complete the second type of communication of the second group of processes in the first column.

5. The heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems according to claim 4 is characterized in that: Except for the first row, the decomposition process of the unit block matrix of each row includes: The process corresponding to the unit block matrix at the diagonal position in the row fills in the unit block matrix at the diagonal position in the row using the basis function information obtained by itself, and updates the filled unit block matrix at the diagonal position in the row using the compression result sent by the second type of communication of the second group of processes and the compression result of all unit block matrices directly above the unit block matrix at the diagonal position in the row; performs LU decomposition on the updated unit block matrix at the diagonal position in the row to obtain the corresponding L matrix and U matrix; and sends the decomposed L matrix to the remaining processes in the first group of processes to complete the first type of communication of the first group of processes in the row; According to the order of the unit block matrices of the row, the process corresponding to the unit block matrix at each non-diagonal position fills in the unit block matrices at all non-diagonal positions corresponding to the row using the basis function information obtained by itself, and updates the unit block matrices at all non-diagonal positions corresponding to the row filled in by the CPU end on the GPU end using the compression results of the unit block matrices at all non-diagonal positions corresponding to the row filled in and the unit block matrices directly above each, as well as the compression results sent by the second type of communication of the second group of processes before; decomposes the unit block matrices at the non-diagonal positions updated by itself using the L matrix obtained by the first type of communication of the first group of processes in the row to obtain the corresponding U matrix, and compresses the obtained U matrix; and, if the unit block matrix at the non-diagonal position of the row is located directly above a unit block matrix at a diagonal position in the block matrix Z, after compressing the obtained U matrix, the current compression result and the compression results of all the unit block matrices directly above are also sent to the matching process in the second group of processes to complete the second type of communication of the first group of processes in the row.

6. The heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems according to claim 4, characterized in that: The decomposition process of the unit block matrix of each column except the first column includes: The process corresponding to the unit block matrix at the diagonal position in the column fills in the unit block matrix at the diagonal position in the column using the basis function information obtained by itself, and updates the filled unit block matrix at the diagonal position in the column using the compression result sent by the second type of communication of the first group of processes and the compression result of all unit block matrices to the left of the unit block matrix at the diagonal position in the column; performs LU decomposition on the updated unit block matrix at the diagonal position in the column to obtain the corresponding L matrix and U matrix; and sends the decomposed U matrix to the remaining processes in the second group of processes to complete the first type of communication of the second group of processes in the column; According to the order of the unit block matrices of the column, the process corresponding to the unit block matrix at each non-diagonal position fills in the unit block matrices at all non-diagonal positions corresponding to itself in the column using the basis function information obtained by itself, and updates the unit block matrices at all non-diagonal positions corresponding to itself in the row filled in by the CPU end on the GPU end using the compression results of the unit block matrices at all non-diagonal positions corresponding to itself in the filled column and the unit block matrices to the immediate left of each, as well as the compression results sent by the second type of communication of the first group of processes before; decomposes the unit block matrices at the non-diagonal positions updated by itself using the U matrix obtained by the first type of communication of the second group of processes in the column to obtain the corresponding L matrix, and compresses the obtained L matrix; and, if the unit block matrix at the non-diagonal position of the column is located to the immediate left of a unit block matrix at a diagonal position in the block matrix Z, after compressing the obtained L matrix, the current compression result and the compression results of all unit block matrices to the immediate left are also sent to the matching process in the first group of processes to complete the second type of communication of the second group of processes in the column.

7. The heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems according to claim 5, characterized in that: For each row other than the first row, the process corresponding to each unit block matrix at each non-diagonal position fills in the unit block matrices at all non-diagonal positions corresponding to the row using the basis function information obtained by itself, and updates the unit block matrices at all non-diagonal positions corresponding to the row filled in by the CPU end on the GPU end using the compression results of the unit block matrices at all non-diagonal positions corresponding to the row filled in by itself and the unit block matrices directly above each other, as well as the compression results sent by the second type of communication of the second group of processes, including: For the process corresponding to the unit block matrix at the non-diagonal position in the row, the filled unit block matrices at all non-diagonal positions corresponding to the row are organized into Local_C data corresponding to the process by columns, the compression results of the unit block matrices directly above the unit block matrices at all non-diagonal positions corresponding to the row are organized into Local_B data corresponding to the process by columns, the compression results of the second type of communication sent by the second group of processes are organized into Local_A data corresponding to the process by columns, the Local_C data, Local_B data and Local_A data corresponding to the process are uploaded to the GPU end, the updated Local_C data is calculated on the GPU end based on the three-step GPU heterogeneous multiplication and downloaded to the CPU end, and the updated results of the unit block matrices at all non-diagonal positions corresponding to the process in the row are obtained after splitting by columns.

8. The heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems according to claim 6, characterized in that: For each column other than the first column, the process corresponding to each unit block matrix at each non-diagonal position fills in all the unit block matrices at non-diagonal positions corresponding to itself in the column using the basis function information obtained by itself, and uses the compression results of all the unit block matrices at non-diagonal positions corresponding to itself in the filled column and the unit block matrices to the left of each other, as well as the compression results sent by the second type of communication of the first group of processes, to update the unit block matrices at all the non-diagonal positions corresponding to itself in the row filled in by the CPU end on the GPU end, including: For the process corresponding to the unit block matrix at the non-diagonal position in the column, the filled unit block matrices at all non-diagonal positions corresponding to the column are organized into rows to form Local_C data corresponding to the process, and the compression results of the unit block matrices to the left of all the unit block matrices at the non-diagonal positions corresponding to the column are organized into rows to form Local_B data corresponding to the process, and the compression results of the second type of communication sent by the second group of processes are organized into rows to form Local_A data corresponding to the process, and the Local_C data, Local_B data and Local_A data corresponding to the process are uploaded to the GPU end, and the updated Local_C data is calculated on the GPU end based on the three-step GPU heterogeneous multiplication and downloaded to the CPU end, and the updated results of the unit block matrices at all non-diagonal positions corresponding to the process in the column are obtained after splitting by rows.

9. The heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems according to claim 7 or 8, characterized in that: In the three-step GPU heterogeneous multiplication, the first two steps are calculated using the cublaszgemm function provided by the cuBLAS platform, and the third step uses a pre-written kernel function for heterogeneous computing.

10. The heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems according to claim 1, characterized in that: The CPU side solves the current coefficient matrix by back-substitution of the matrix equations according to the decomposition result of the block matrix Z, including: The CPU obtains the L matrix and the U matrix decomposed from the block matrix Z from the decomposition result of the block matrix Z, and calculates the excitation matrix V according to the formula LUI=V based on the original solution formula ZI=V of the moment method; wherein I represents the current coefficient matrix to be solved; UI is set as the intermediate process matrix Y, and the intermediate process matrix Y is calculated according to the formula LY=V; The current coefficient matrix I is calculated according to the formula UI=Y.

Images