SIMD (Single Instruction Multiple Data)-oriented data processing method for parallel iterative solution

Abstract

The invention provides a data processing method for SIMD (Single Instruction Multiple Data)-oriented parallel iterative solution, which comprises the following steps: acquiring to-be-processed data (including 2m pattern matrixes in a complex form), storing each pattern matrix in a real number newCSR format (including a ValuesR array, a ValuesI array, a rowPtr array and a colIndices array), and vectorizing to obtain a vectorized fused newCSR format (including a newValuesR array, a newValuesI array, a rowPtr array and a colIndices array), and reading vectorized new CSR format fused data by using an SIMD extension component, carrying out parallel iterative solution, and outputting an iterative solution result. In this way, parallel iterative solution calculation can be carried out through the SIMD extension component in the SIMD architecture, and accelerated solution is achieved.

Description

Technical Field

[0001] The present application relates to the field of data processing technology, and in particular to a data processing method for SIMD-oriented parallel iterative solution. Background Art

[0002] In the field of numerical analysis (such as electromagnetic field numerical analysis), there is often a need to solve matrices. Traditional solutions include direct solution method, iterative solution method, etc.

[0003] Iterative methods for linear equations are a class of methods used to solve large sparse linear equations by gradually approximating the exact solution. These methods are usually superior to direct methods in terms of computational efficiency and memory usage, especially when dealing with large-scale problems. The following are several common iterative methods for linear equations and their brief introduction:

[0004] 1. Jacobi Iteration

[0005] Basic idea: decompose the linear equation system Ax = b into Dx + Lx + Ux = b, where D is a diagonal matrix, L is a lower triangular matrix, and U is an upper triangular matrix. Then, by iterating the formula x (k+1) =D -1 (b-(L+U)x (k) ) to gradually approach the solution.

[0006] Features: Simple and easy to implement, but slow convergence speed, suitable for diagonally dominant matrices.

[0007] 2. Gauss-Seidel Iteration

[0008] Basic idea: Similar to Jacobi iteration, but the latest component value is used immediately in each iteration. The iteration formula is x (K+1) =(D+L) -1 (b-Ux (k) ).

[0009] Features: Usually converges faster than the Jacobi iteration method, but the computational complexity is slightly higher.

[0010] 3. Successive Over-Relaxation (SOR)

[0011] Basic idea: Based on the Gauss-Seidel iteration method, a relaxation parameter ω is introduced to accelerate convergence by adjusting ω. The iteration formula is:

[0012]

[0013] Features: By choosing a suitable ω, the convergence speed can be significantly accelerated, but improper selection may lead to divergence.

[0014] 4. Conjugate Gradient (CG)

[0015] Basic idea: Applicable to symmetric positive definite matrix A, the solution is gradually approached by constructing a set of conjugate directions. Only matrix-vector multiplication is required during the iteration, and there is no need to store the entire matrix.

[0016] Features: fast convergence speed, high computational efficiency, especially suitable for large sparse matrices.

[0017] 5. Biconjugate Gradient (BCG)

[0018] Basic idea: Applicable to asymmetric matrices, and handles asymmetric situations by constructing two mutually conjugated sequences.

[0019] Features: Works well for nonsymmetric matrices, but may encounter convergence problems, especially if the matrix is ​​ill-conditioned or has complex eigenvalues.

[0020] 6. Biconjugate Gradient Stabilized (BiCGSTAB)

[0021] Basic idea: An additional stabilization step is introduced on the basis of the BCG method to improve convergence and stability.

[0022] Features: More stable than BCG, converges faster, and is suitable for asymmetric matrices.

[0023] 7. Generalized Minimal Residual (GMRES)

[0024] Basic idea: gradually approximate the solution by minimizing the residual norm, which is applicable to non-symmetric matrices. Each step reconstructs a Krylov subspace.

[0025] Features: Good convergence, but high computational complexity and the need to store multiple vectors.

[0026] 8. Preconditioned Conjugate Gradient (PCG)

[0027] Basic idea: Based on the conjugate gradient method, a preconditioner M is introduced. -1 A to accelerate convergence.

[0028] Features: By selecting a suitable preconditioner, the convergence speed can be significantly improved. It is applicable to symmetric positive definite matrices.

[0029] These methods have their own advantages and disadvantages. Choosing the right method depends on the specific properties of the problem, such as the symmetry, sparsity, and degree of pathology of the matrix.

[0030] Modern processor architectures (mainstream ones include x86 architecture and ARM architecture, ARM NEON is the SIMD instruction set under the ARM platform, the x86 platform has MMX, SSE, AVX and other instruction sets, and some other architectures also have some SIMD instruction sets) can support the simultaneous execution of multiple instructions through their SIMD extension components (vectorization components). SIMD, which stands for single instruction multiple data, means that a single operation instruction can execute multiple data streams, thereby increasing the operation speed. SIMD extension components refer to dedicated hardware units in the processor for performing SIMD operations. These extension components usually include wide registers, dedicated SIMD instruction sets and corresponding execution units. The following are several common SIMD extension components and their characteristics.

[0031] 1. Intel's SIMD extensions: name SSE (Streaming SIMD Extensions); versions: SSE, SSE2, SSE3, SSSE3, SSE4.1, SSE4.2; register width: 128 bits; data types: single-precision floating point (32 bits), double-precision floating point (64 bits), integers (8 bits, 16 bits, 32 bits, 64 bits); uses: widely used in floating-point operations, multimedia processing, and scientific computing; registers: XMM0 to XMM15 (possibly more on newer processors).

[0032] 2. AMD's SIMD extension: Name: 3DNow!; Versions: 3DNow!, 3DNow! Professional, 3DNow! +, Enhanced 3DNow!; Register width: 64 bits; Data type: single-precision floating point (32 bits); Purpose: mainly used for floating-point operations, especially for games and multimedia applications; Registers: MM0 to MM7.

[0033] 3. ARM's SIMD extension: Name: NEON register; Width: 64 bits and 128 bits; Data types: single-precision floating point (32 bits), double-precision floating point (64 bits), integer (8 bits, 16 bits, 32 bits, 64 bits); Uses: widely used in multimedia processing, signal processing and high-performance computing; Registers: Q0 to Q31 (128 bits), D0 to D31 (64 bits).

[0034] Take NEON as an example. NEON technology is an advanced SIMD (single instruction multiple data) architecture for ArmCortex-A series processors. It can accelerate multimedia and signal processing algorithms such as video encoders / decoders, 2D / 3D graphics, games, audio and voice processing, image processing, phone calls and sounds. NEON instructions perform "packed SIMD" processing, registers are considered vectors of elements of the same data type, and data types support: signed / unsigned 8-bit, 16-bit, 32-bit, 64-bit.

[0035] There are 32 NEON vector registers in ARMv8, with a length of 128 bits. They can be used to process both scalar operands and vector operands. These registers can store multiple operands of the same data type, such as two 64-bit double type operands or four 32-bit float type operands. During the execution of NEON instructions, multiple operands on the same vector register are processed simultaneously, thereby achieving parallel processing of different data. Using vector registers to accelerate programs requires reading data from memory into them in advance, so the continuity of the accessed data in memory has an important impact on the performance of vector registers, and data alignment can significantly improve the performance of vector registers.

[0036] In some fields represented by electromagnetic computing, we are often faced with solving equations with the same distribution of non-zero elements. For example, in the field of electromagnetic computing, when the model and mesh division remain unchanged, the distribution of non-zero elements in the equations generated at different frequencies is the same. Since there are a large number of complex matrices (dense complex matrices, sparse complex matrices, etc.) in the field of electromagnetic numerical analysis that need to be solved linearly, the current technology requires calling a solution function for each solution, and there is no way to solve two or more linear equations in parallel.

[0037] Based on this, the inventors of the present application have designed a parallel iterative solution for SIMD, which can realize parallel iterative solution of pattern matrices (large complex matrices with the same non-zero element distribution pattern) in various scenarios to achieve accelerated solution. Summary of the invention

[0038] The purpose of the embodiments of the present application is to provide a data processing method for parallel iterative solution of SIMD, so as to perform parallel iterative solution calculation through SIMD extension components to achieve accelerated solution.

[0039] In order to achieve the above purpose, the embodiments of the present application are implemented in the following manner:

[0040] In a first aspect, an embodiment of the present application provides a data processing method for SIMD parallel iterative solution, comprising: obtaining data to be processed, wherein the data to be processed includes 2 m pattern matrix, m∈Z + , each pattern matrix has the same non-zero element distribution pattern and is a complex matrix; each pattern matrix is ​​stored in the real number new_CSR format, where the real number new_CSR format contains the ValuesR array, the ValuesI array, the rowPtr array, and the colIndices array. The ValuesR array is used to store the real part of the non-zero elements of the pattern matrix, the ValuesI array is used to store the imaginary part of the non-zero elements of the pattern matrix, the rowPtr array is used to store the position of the first non-zero element of each row in the pattern matrix in the values ​​array, and the colIndices array is used to store the column index of each non-zero element in the pattern matrix; for 2 m The real number new_CSR format of the pattern matrix is ​​vectorized to obtain the vectorized fused new_CSR format, where the fused new_CSR format includes the new_ValuesR array, the new_ValuesI array, the rowPtr array and the colIndices array, and the new_ValuesR array includes several real part arrays, each of which is used to store 2 m The real part of the non-zero elements at the same position of the pattern matrix, the new_ValuesI array contains the corresponding number of imaginary arrays, each imaginary array is used to store 2 m The imaginary part of the non-zero elements of the pattern matrices at the same position is obtained; the SIMD extension component is used to read the vectorized fused new_CSR format data, perform parallel iterative solutions, and output the iterative solution results.

[0041] In combination with the first aspect, in a first possible implementation method of the first aspect, each pattern matrix is ​​stored in a real number new_CSR format, including: storing each pattern matrix in a complex number CSR format, wherein the complex number CSR format includes a Values ​​array, a rowPtr array, and a colIndices array, the Values ​​array is used to store the non-zero elements of the pattern matrix, the rowPtr array is used to store the position of the first non-zero element of each row in the pattern matrix in the values ​​array, and the colIndices array is used to store the column index of each non-zero element in the pattern matrix; for the complex number CSR format of each pattern matrix: splitting the complex number non-zero elements of the Values ​​array in the complex number CSR format into the real part of the non-zero elements and the imaginary part of the non-zero elements, storing the real parts of the non-zero elements in sequence to obtain the ValuesR array, storing the imaginary parts of the non-zero elements in sequence to obtain the ValuesI array, and the rowPtr array and the colIndices array remain unchanged.

[0042] In combination with the first aspect, in a second possible implementation manner of the first aspect, m The real number new_CSR format of the pattern matrix is ​​vectorized to obtain the vectorized fusion new_CSR format, including: m The real parts of the non-zero elements of the ValuesR array at the same position in the new_CSR format of the pattern matrix are combined into a 2-bit array in the order of the pattern matrix. m The real part array of the real part of the non-zero elements is obtained, and the new_ValuesR array containing several real part arrays is obtained; 2 m The imaginary parts of the non-zero elements in the ValuesI array at the same position in the real number new_CSR format of the pattern matrix are combined into a 2 m The imaginary part array of the non-zero elements of the imaginary part is obtained, and the new_ValuesR array containing several imaginary part arrays is obtained; the rowPtr array remains unchanged; the colIndices array remains unchanged.

[0043] In combination with the first aspect, in a third possible implementation of the first aspect, using a SIMD extension component to read vectorized fused new_CSR format data and perform parallel iterative solving includes: using a SIMD extension component to read 2 m The fused new_CSR format quaternion of the pattern matrix, the fused new_CSR format right value Breal and Bimag, the dimension n of the pattern matrix, the maximum number of iterations IMAX, and the convergence tolerance TOL, where the fused new_CSR format quaternion is 2 mThe matrix A in the pattern matrix Ax=b is stored in the real number new_CSR format and then vectorized, including the new_ValuesR array, new_ValuesI array, rowPtr array and colIndices array, and the right value of the fused new_CSR format is 2 m The right-hand vector b in the pattern matrix Ax=b is stored in the real number new_CSR format and then vectorized, including Breal and Bimag, Breal is the real part of the right value, Bimag is the imaginary part of the right value; initialize the intermediate vectors AXR and AXI, search direction vectors Preal and Pimag, residual vectors Rreal and Rimag, solution vectors XR and XI, residual norm VNRM, errors RSS1,…,RSSi…,RSS2 m , and assign values: Rreal=Breal,Rimag=Bima,Preal=Rreal,Pimag=Pimag,RSS1=…RSSi=…RSS2 m =1.0,i∈[1,2 m ], where AXR and AXI are the real and imaginary parts of the intermediate vector, Preal and Pimag are the real and imaginary parts of the search direction vector, Rreal and Rimag are the real and imaginary parts of the residual vector, and RSSi is the error corresponding to the i-th pattern matrix; based on the intermediate vectors AXR and AXI, the search direction vectors Preal and Pimag, the residual vectors Rreal and Rimag, the residual norm VNRM, the errors RSS1,…,RSSi…,RSS2 m , use the SIMD extension component to iteratively update the solution vectors XR and XI, and finally output the iterative solution vectors XR and XI; decompose the solution vectors XR and XI to obtain 2 m The solution vector of the pattern matrix.

[0044] In combination with the third possible implementation of the first aspect, in a fourth possible implementation of the first aspect, based on the intermediate vectors AXR and AXI, the search direction vectors Preal and Pimag, the residual vectors Rreal and Rimag, the residual norm VNRM, the errors RSS1, ..., RSSi ..., RSS2 m, use the SIMD extension component to iteratively update the solution vectors XR and XI, and finally output the iterative solution vectors XR and XI, including: S1: determine whether the current iteration termination condition is met; S2: if the iteration termination condition is met, end the iteration and output the iterative solution vectors XR and XI; S3: if the iteration termination condition is not met, start this iteration: S31: call the improved cnorm function to record the residual vector inner product value RNORM under the current iteration number; S32: call the improved productAll function to calculate the intermediate vectors AXR and AXI under the current iteration number; S33: based on the residual vector inner product value RNORM and the intermediate vectors AXR and AXI under the current iteration number, calculate the update step size α under the current iteration number; S34: based on the update step size α under the current iteration number and the search direction vector P real and Pimag, update the solution vectors XR and XI at the current number of iterations; S35: based on the update step size α and the intermediate vectors AXR and AXI at the current number of iterations, update the residual vectors Rreal and Rimag, and obtain the residual vectors Rreal′ and Rimag′ at the current number of iterations; S36: based on the residual vector inner product value RNORM at the current number of iterations, the residual vectors Rreal′ and Rimag′ at the current number of iterations, calculate the conjugate parameter γ; S37: based on the conjugate parameter γ, update the search direction vectors Preal and Pimag; S38: add 1 to the number of iterations, calculate the residual norm VNRM2 after the iteration, and calculate the errors RSS1,…,RSSi…,RSS2 at the current number of iterations based on the initialized residual norm VNRM and the residual norm VNRM2 after the iteration. m , and then return to S1.

[0045] In combination with the fourth possible implementation of the first aspect, in a fifth possible implementation of the first aspect, the improved cnorm function is called to record the residual vector inner product value RNORM at the current number of iterations, including: initializing vectors sum1 and sum2; the SIMD extension component reads 2 in the residual vector Rreal; m data, read 2 of the residual vector Rimag m data, where the SIMD extension unit reads data in steps of 2 m , Rreal is the real part of the residual vector R before updating at the current number of iterations, and Rimag is the imaginary part of the residual vector R before updating at the current number of iterations; Calculate according to the following formula2 m The residual vector inner product value RNORM:

[0046] RNORM=R T R,

[0047] RNORMreal=RrealT real-Rimag T Rimag,

[0048] RNORMimag=Rreal T Rimgag+Rimag T Rreal,

[0049] Among them, RNORMreal is the real part of the residual norm RNORM, RNORMimag is the imaginary part of the residual vector inner product value RNORM, Rreal T is the transpose of Rreal, Rimag T is the transpose of Rimag; by CCR[0]=sum1[0], CCR[1]=sum1[1], ..., CCR[k]=sum1[k], ..., CCR[2 m -1]=sum1[2 m -1] Save 2 m The real part of the inner product value of the residual vector RNORMreal is obtained by CCR[[2 m ]=sum2[0]、CCR[[2 m +1] = sum2[1], ..., CCR[2 m + k] = sum2[k], ..., CCR[2 m+1 -1]=sum2[2 m -1] Save 2 m The imaginary part of the inner product value of the residual vector RNORMimag.

[0050] In combination with the fourth possible implementation of the first aspect, in a sixth possible implementation of the first aspect, the improved productAll function is called to calculate the intermediate vectors AXR and AXI under the current number of iterations, including: the SIMD extension component reads the search direction vectors Preal and Pimag, and fuses the quadruple in the new_CSR format: new_ValuesR array, new_ValuesI array, rowPtr array and colIndices array, wherein the step size of the SIMD extension component reading data is 2 m ; Calculate the intermediate vectors AXR and AXI at the current iteration number according to the following formula:

[0051] AXR=AR·Preal-AI*Pimag,

[0052] AXI=AR*Pimag+AI*Preal,

[0053] Among them, AXR is the real part of the intermediate vector AX, AR is the real part of the pattern matrix, which is determined by reading the new_ValuesR array, rowPtr array and colIndices array through the SIMD extension component, Preal is the real part of the search direction vector P, AXI is the imaginary part of the intermediate vector AI, AI is the imaginary part of the pattern matrix, which is determined by reading the new_ValuesI array, rowPtr array and colIndices array through the SIMD extension component, and Pimag is the imaginary part of the search direction vector P.

[0054] In combination with the fifth possible implementation of the first aspect, in the seventh possible implementation of the first aspect, based on the residual vector inner product value RNORM and the intermediate vectors AXR and AXI at the current number of iterations, the update step size α at the current number of iterations is calculated, including: the SIMD extension component reads the residual vector inner product value RNORMreal and RNORMimag at the current number of iterations, the search direction vectors Preal and Pimag, and the intermediate vectors AXR and AXI at the current number of iterations, wherein the step size of the SIMD extension component reading data is 2 m ;

[0055] The update step size α at the current number of iterations is calculated according to the following formula:

[0056]

[0057] Further:

[0058]

[0059] Among them, αr is the real part of the update step size α, and αi is the imaginary part of the update step size α.

[0060] In combination with the seventh possible implementation of the first aspect, in an eighth possible implementation of the first aspect, based on the update step size α and the search direction vectors Preal and Pimag at the current iteration number, the solution vectors XR and XI at the current iteration number are updated, including: the SIMD extension component reads the update step size α and the search direction vectors Preal and Pimag at the current iteration number, wherein the step size of the SIMD extension component reading data is 2 m ; Calculate and update the solution vector at the current iteration number according to the following formula:

[0061] XR′=XR+ar*Preal-αi*Pimag,

[0062] XI′=XI+αr*Pimag+ai*Preal,

[0063] Among them, XR′ is the real part of the solution vector X after update under the current iteration number, XR is the real part of the solution vector X before update under the current iteration number, that is, the real part of the solution vector X updated in the last iteration, αr is the real part of the update step α, Preal is the real part of the search direction vector P, XI′ is the imaginary part of the solution vector X after update under the current iteration number, XI is the imaginary part of the solution vector X before update under the current iteration number, that is, the imaginary part of the solution vector X updated in the last iteration, αi is the imaginary part of the update step α, and Pimag is the imaginary part of the search direction vector P.

[0064] In combination with the seventh possible implementation of the first aspect, in a ninth possible implementation of the first aspect, based on the update step size α and the intermediate vectors AXR and AXI at the current number of iterations, the residual vectors Rreal and Rimag are updated to obtain the residual vectors Rreal′ and Rimag′ at the current number of iterations, including: a SIMD extension component reads the update step size α and the intermediate vectors AXR and AXI at the current number of iterations, wherein the step size of the SIMD extension component reading the data is 2 m ; Update the residual vector at the current iteration number according to the following formula:

[0065] Rreal′=Rreal-αr*AXR+αi*AXI,

[0066] Rimag′=Rimag-αr*AXI-αi*AXR,

[0067] Among them, Rreal′ is the real part of the residual vector R after update at the current iteration number, Rreal is the real part of the residual vector R before update at the current iteration number, that is, the real part of the residual vector R updated in the previous iteration, αr is the real part of the update step α, AXR is the real part of the intermediate vector AX at the current iteration number, Rimag′ is the imaginary part of the residual vector R after update at the current iteration number, Rimag is the imaginary part of the residual vector R before update at the current iteration number, that is, the imaginary part of the residual vector R updated in the previous iteration, αi is the imaginary part of the update step α, and AXI is the imaginary part of the intermediate vector AX at the current iteration number.

[0068] In combination with the ninth possible implementation of the first aspect, in the tenth possible implementation of the first aspect, based on the residual vector inner product value RNORM at the current number of iterations, the residual vectors Rreal′ and Rimag′ at the current number of iterations, the conjugate parameter γ is calculated, including: a SIMD extension component reads the residual vector inner product value RNORM at the current number of iterations, the residual vectors Rreal′ and Rimag′ at the current number of iterations, wherein the step size of the SIMD extension component reading the data is 2 m ; Calculate the conjugate parameter γ according to the following formula:

[0069]

[0070] Among them, γ is the conjugate parameter, which includes the real part γr and the imaginary part γi, as follows:

[0071]

[0072] Among them, Rreal′ is the real part of the residual vector R′ after updating at the current number of iterations, Rreal′ T is the transpose of Rreal′, Rimag′ is the imaginary part of the residual vector R′ updated at the current number of iterations, and Rimag′ T is the transpose of Rimag′, Rreal is the real part of the residual vector R before updating at the current iteration number, that is, the real part of the residual vector R updated in the last iteration, Rreal T is the transpose of Rreal, Rimag is the imaginary part of the residual vector R before update at the current iteration number, that is, the imaginary part of the residual vector R updated in the last iteration, Rimag T is the transpose of Rimag.

[0073] In combination with the tenth possible implementation of the first aspect, in the eleventh possible implementation of the first aspect, the search direction vectors Preal and Pimag are updated based on the conjugate parameter γ, including: the SIMD extension component reads the search direction vectors Preal and Pimag, the real part γr and the imaginary part γi of the conjugate parameter γ at the current number of iterations, and the residual vectors Real′ and Rimag′ at the current number of iterations, wherein the step length of the SIMD extension component reading the data is 2 m ; Calculate and update the search direction vectors Preal and Pimag according to the following formula:

[0074] Preal′=Rreal′+γr*Preal-γi*Pimag,

[0075] Pimag′=Rimag′+γr*Pimag+γl*Preal,

[0076] Among them, Preal′ is the real part of the search direction vector P′ after update at the current iteration number, Preal is the real part of the search direction vector P before update at the current iteration number, Pimag′ is the imaginary part of the search direction vector P′ after update at the current iteration number, and Pimag is the imaginary part of the search direction vector P before update at the current iteration number.

[0077] In combination with the ninth possible implementation manner of the first aspect, in a twelfth possible implementation manner of the first aspect, the residual norm VNRM2 after iteration is calculated, and based on the initialized residual norm VNRM and the residual norm VNRM2 after iteration, the error RSS1,…,RSSi…,RSS2 at the current number of iterations is calculated. m , including: the SIMD extension component reads the residual vectors Rreal′ and Rimag′ at the current number of iterations, where the step length of the SIMD extension component reading data is 2 m ; Calculate the residual norm VNRM2 of each pattern matrix after this iteration according to the following formula:

[0078]

[0079] Where VNRM2 is the residual norm calculated by the pattern matrix at the current number of iterations, R′ is the updated residual vector at the current number of iterations, n is the dimension of the number of rows or columns in the pattern matrix, and Rreal′ j is the jth value of Rreal′, Rimag′ j is the kth value of Rimag′; the SIMD extension component reads the initialized residual norm VNRM and the iterated residual norm VNRM2, where the step length of the SIMD extension component reading data is 2 m ; Calculate the error of each pattern matrix at the current number of iterations according to the following formula:

[0080]

[0081] Among them, RSSi is the error of the i-th pattern matrix at the current number of iterations, i∈[1,2 m ].

[0082] In combination with the twelfth possible implementation method of the first aspect, in the thirteenth possible implementation method of the first aspect, the iteration termination condition is: the number of iterations iter reaches the maximum number of iterations IMAX, or the error RSS of each pattern matrix at the current number of iterations is not greater than the convergence tolerance TOL.

[0083] Beneficial effects:

[0084] 1. In the field of numerical analysis (such as electromagnetic numerical analysis, mechanical numerical analysis, etc.), there are a large number of complex sparse matrices that need to be solved linearly, but each solution requires calling a solution function, and there is no way to solve two or more linear equations in parallel. This solution uses SIMD extension components to design 2 m A parallel iterative solution for the pattern matrix is ​​proposed to achieve accelerated solution. Specifically, the data to be processed (including 2 mpattern matrix, m∈Z + , each pattern matrix has the same non-zero element distribution pattern and is a complex matrix); store each pattern matrix in the real number new_CSR format, where the real number new_CSR format contains the ValuesR array, the ValuesI array, the rowPtr array, and the colIndices array. The ValuesR array is used to store the real part of the non-zero elements of the pattern matrix, the ValuesI array is used to store the imaginary part of the non-zero elements of the pattern matrix, the rowPtr array is used to store the position of the first non-zero element of each row in the pattern matrix in the values ​​array, and the colIndices array is used to store the column index of each non-zero element in the pattern matrix; for 2 m The real number new_CSR format of the pattern matrix is ​​vectorized to obtain the vectorized fused new_CSR format, where the fused new_CSR format includes the new_ValuesR array, the new_ValuesI array, the rowPtr array and the colIndices array, and the new_ValuesR array includes several real part arrays, each of which is used to store 2 m The real part of the non-zero elements at the same position of the pattern matrix, the new_ValuesI array contains the corresponding number of imaginary arrays, each imaginary array is used to store 2 m The imaginary part of the non-zero elements of the same position of the pattern matrix; use the SIMD extension component to read the vectorized fusion new_CSR format data, perform parallel iterative solution, and output the iterative solution result. This solution improves the storage format of CSR and vectorizes the corresponding data of multiple pattern matrices so that the SIMD extension component (such as NEON register) can read data and execute SIMD instructions (such as NEON instructions) in parallel, improves the relevant iterative solution algorithm, and thus achieves 2 m The parallel iterative solution of the pattern matrix greatly improves the iterative solution efficiency.

[0085] 2. This solution can be achieved by storing each pattern matrix in a complex CSR format. The complex CSR format contains a Values ​​array, a rowPtr array, and a colIndices array. The Values ​​array is used to store the non-zero elements of the pattern matrix, the rowPtr array is used to store the position of the first non-zero element of each row in the pattern matrix in the values ​​array, and the colIndices array is used to store the column index of each non-zero element in the pattern matrix; for the complex CSR format of each pattern matrix: split the complex non-zero elements of the Values ​​array in the complex CSR format into the real part of the non-zero elements and the imaginary part of the non-zero elements, store the real part of the non-zero elements in sequence to obtain the ValuesR array, store the imaginary part of the non-zero elements in sequence to obtain the ValuesI array, and the rowPtr array and the colIndices array remain unchanged. In this way, the existing complex CSR format can be used for further splitting and processing to obtain a new real number new_CSR format, which can make full use of existing technologies and be applicable to a wide range of application scenarios. And 2 m The real parts of the non-zero elements of the ValuesR array at the same position in the new_CSR format of the pattern matrix are combined into a 2-bit array in the order of the pattern matrix. m The real part array of the real part of the non-zero elements is obtained, and the new_ValuesR array containing several real part arrays is obtained; 2 m The imaginary parts of the non-zero elements in the ValuesI array at the same position in the real number new_CSR format of the pattern matrix are combined into a 2 m The imaginary part array of the imaginary part of the non-zero elements is obtained to obtain the new_ValuesR array containing several imaginary part arrays; the rowPtr array remains unchanged; the colIndices array remains unchanged. Since the use of vector registers to accelerate the program requires the data in the memory to be read into it in advance, vectorization is performed based on the real number new_CSR format to obtain the fused new_CSR format, which can achieve vectorization and make the fused new_CSR format adapt to the SIMD extension component. Data alignment can significantly improve the performance of the SIMD extension component, realize parallel iterative solution calculation, effectively improve the calculation speed of parallel processing, greatly improve calculation efficiency, and shorten iterative solution time.

[0086] 3. The parallel iterative solution of this solution improves the iterative solution algorithm so that the entire iterative solution algorithm can be applied to 2 mThe parallel iterative solution calculation of the pattern matrix is ​​designed, and the corresponding improved formula is designed to realize parallel iterative update. At the same time, the operation step size of the SIMD extension component (such as reading data) is adjusted to ensure the stable operation of the parallel iterative solution. The test case uses this solution (taking the parallel accelerated iterative solution of NEON register as an example) to iteratively solve two complex sparse matrices of size 122036*122036, where the non-zero value of each matrix is ​​8168600. The comparison method is not to use the SIMD extension component to accelerate the iterative solution, and perform CG iterative solution on the two matrices respectively. The solution time is 2696.53 seconds and 8087.68 seconds respectively, and the iterative solution efficiency is improved by 66.66%. If more matrices are solved, the acceleration performance will be more obvious.

[0087] In order to make the above-mentioned objects, features and advantages of the present application more obvious and easy to understand, preferred embodiments are specifically cited below and described in detail with reference to the attached drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0088] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings required for use in the embodiments of the present application will be briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present application and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without paying creative work.

[0089] Figure 1 A flowchart of a data processing method for SIMD parallel iterative solution provided in an embodiment of the present application.

[0090] Figure 2 A schematic diagram of storing the pattern matrix in the real number new_CSR format.

[0091] Figure 3 A schematic diagram of storing the pattern matrix in complex CSR format.

[0092] Figure 4 A schematic diagram of splitting the Values ​​array in the complex CSR format into the ValuesR array and the ValuesI array.

[0093] Figure 5 A schematic diagram of vectorizing the real number new_CSR format of two pattern matrices to obtain a fused new_CSR format.

[0094] Figure 6 The overall flow chart of the improved iterative solution algorithm. DETAILED DESCRIPTION

[0095] The technical solutions in the embodiments of the present application will be described below in conjunction with the drawings in the embodiments of the present application.

[0096] In the field of numerical analysis (such as electromagnetic field numerical analysis, mechanical numerical analysis, etc.), there is often a need to solve matrices. Among these matrices that need to be solved (such as large dense complex matrices, large sparse complex matrices, etc.), there are many pattern matrices, where the pattern matrix represents a matrix with the same non-zero element distribution pattern. For example, in the field of electromagnetic calculations, when the model and mesh division remain unchanged, the distribution positions of the non-zero elements of the equations generated at different frequencies are the same. If these matrices are to be solved iteratively, each solution requires a call to the solver function, and there is no way to solve two or more linear equations in parallel.

[0097] This embodiment designs a data processing method for SIMD parallel iterative solution for such large complex matrices with the same non-zero element distribution pattern (actually, it can be used not only for large complex matrices, but also for iterative solution of all such pattern matrices. Generally, small matrices can be solved directly with higher efficiency). By using SIMD extension components, the simultaneous horizontal solution of multiple equation groups with the same non-zero element distribution is realized. For example, a model in the electromagnetic field has been grid discretized, and the electromagnetic responses at two frequencies of 100 Hz and 10 Hz are solved. Each equation group needs to be solved, and the non-zero element distribution positions of these two equation groups are exactly the same, resulting in two pattern matrices. Therefore, it is suitable to use the method proposed in the present invention to horizontally vectorize the solution process to achieve parallel solution.

[0098] See also Figure 1 , Figure 1 Flow chart of the data processing method for SIMD parallel iterative solution provided in this embodiment. In this embodiment, the data processing method for SIMD parallel iterative solution may include step S10, step S20, step S30, and step S40.

[0099] First, step S10 may be executed.

[0100] Step S10: Obtain the data to be processed, wherein the data to be processed includes 2 m pattern matrix, m∈Z + , each pattern matrix has the same non-zero element distribution pattern and is a complex matrix.

[0101] In this embodiment, the data to be processed, that is, the pattern matrix to be solved (usually also including boundary conditions, written in the form of Ax=b), can be obtained. The pattern matrix to be processed once needs to satisfy 2 m The number of +, each pattern matrix has the same non-zero element distribution pattern, and is a complex matrix. At present, the latest ARM9 platform supports the execution of 4 double-precision instructions at the same time, and it may be expanded in the future. Although the technical difficulty is high and the double-precision instructions that can be processed synchronously are limited (currently the level of executing 2-4 double-precision instructions at the same time), there is still a large application potential, and the idea of ​​this solution can theoretically achieve a higher level of parallel processing, so as to greatly improve processing efficiency. This embodiment is introduced by taking 2 pattern matrices as an example for convenience of explanation, which should not be regarded as a limitation of this application.

[0102] After the data to be processed is obtained, step S20 may be executed.

[0103] Step S20: Store each pattern matrix in the real number new_CSR format, wherein the real number new_CSR format includes the ValuesR array, the ValuesI array, the rowPtr array and the colIndices array, the ValuesR array is used to store the real part of the non-zero elements of the pattern matrix, the ValuesI array is used to store the imaginary part of the non-zero elements of the pattern matrix, the rowPtr array is used to store the position of the first non-zero element of each row in the pattern matrix in the values ​​array, and the colIndices array is used to store the column index of each non-zero element in the pattern matrix.

[0104] In this embodiment, each pattern matrix needs to be stored in a real number new_CSR format, which includes a ValuesR array, a ValuesI array, a rowPtr array, and a colIndices array. Figure 2 As shown, the ValuesR array is used to store the real part of the non-zero elements of the pattern matrix, the ValuesI array is used to store the imaginary part of the non-zero elements of the pattern matrix, the rowPtr array is used to store the position of the first non-zero element of each row in the pattern matrix in the values ​​array, and the colIndices array is used to store the column index of each non-zero element in the pattern matrix.

[0105] Since the prior art already has a CSR storage format (for complex matrices, they are stored in a complex CSR format), the complex CSR format can be further processed and converted into a real number new_CSR format.

[0106] like Figure 3As shown, illustratively, each pattern matrix can be stored in a complex CSR format (which can be done with the existing technology), the complex CSR format includes a Values ​​array, a rowPtr array, and a colIndices array, the Values ​​array is used to store the non-zero elements of the pattern matrix, the rowPtr array is used to store the position of the first non-zero element of each row in the pattern matrix in the values ​​array, and the colIndices array is used to store the column index of each non-zero element in the pattern matrix.

[0107] like Figure 4 As shown, for the complex CSR format of each pattern matrix: the non-zero elements in the complex form of the Values ​​array in the complex CSR format can be split into the real part of the non-zero elements and the imaginary part of the non-zero elements, the real part of the non-zero elements are stored in sequence to obtain the ValuesR array, the imaginary part of the non-zero elements are stored in sequence to obtain the ValuesI array, and the rowPtr array and the colIndices array remain unchanged.

[0108] Thus, the real number new_CSR format of the pattern matrix can be obtained. Of course, in addition to this solution, the pattern matrix can also be stored in the real number new_CSR format according to the rules by directly improving the code during storage, which is not limited here.

[0109] Get 2 m After the real number new_CSR format corresponding to the pattern matrix is ​​obtained, step S30 can be executed.

[0110] Step S30: Pair 2 m The real number new_CSR format of the pattern matrix is ​​vectorized to obtain the vectorized fused new_CSR format, where the fused new_CSR format includes the new_ValuesR array, the new_ValuesI array, the rowPtr array and the colIndices array, and the new_ValuesR array includes several real part arrays, each of which is used to store 2 m The real part of the non-zero elements at the same position of the pattern matrix, the new_ValuesI array contains the corresponding number of imaginary arrays, each imaginary array is used to store 2 m The imaginary parts of the non-zero elements of the pattern matrices at the same position.

[0111] In this embodiment, 2 m The real number new_CSR format of the pattern matrix is ​​vectorized. Specifically, 2 m The real parts of the non-zero elements of the ValuesR array at the same position in the new_CSR format of the pattern matrix are combined into a 2-bit array in the order of the pattern matrix. mThe real part array of the real part of the non-zero elements is obtained, and the new_ValuesR array containing several real part arrays is obtained; 2 m The imaginary parts of the non-zero elements in the ValuesI array at the same position in the real number new_CSR format of the pattern matrix are combined into a 2 m The imaginary part array of the non-zero elements of the imaginary part is obtained, and the new_ValuesR array containing several imaginary part arrays is obtained; the rowPtr array remains unchanged; the colIndices array remains unchanged.

[0112] like Figure 5 As shown, taking the vectorization of the real number new_CSR format corresponding to two pattern matrices as an example, the two pattern matrices are converted into the corresponding real number new_CSR format, and then vectorized to obtain the corresponding fused new_CSR format.

[0113] Get 2 m After the new_CSR format corresponding to the pattern matrices is merged, step S40 can be executed.

[0114] Step S40: Use the SIMD extension component to read the vectorized fused new_CSR format data, perform parallel iterative solution, and output the iterative solution result.

[0115] In order to achieve 2 m The parallel iterative solution of the pattern matrix is ​​designed in this embodiment. The overall process is as follows: Figure 6 shown.

[0116] First, initialize:

[0117] Using SIMD extensions to read 2 m The fused new_CSR format quaternion of the pattern matrix, the fused new_CSR format right value Breal and Bimag, the dimension n of the pattern matrix, the maximum number of iterations IMAX, and the convergence tolerance TOL, where the fused new_CSR format quaternion is 2 m The matrix A in the pattern matrix Ax=b is stored in the real number new_CSR format and then vectorized, including the new_ValuesR array, new_ValuesI array, rowPtr array and colIndices array, and the right value of the fused new_CSR format is 2 m The right-hand side vector b in the pattern matrix Ax=b is obtained by vectorization after being stored in the real number new_CSR format, and includes Breal and Bimag, Breal is the real part of the right value, and Bimag is the imaginary part of the right value.

[0118] Then, initialize the intermediate vectors AXR and AXI, search direction vectors Preal and Pimag, residual vectors Rreal and Rimag, solution vectors XR and XI, residual norm VNRM, errors RSS1,…,RSSi…,RSS2 m , and assign values: Rreal=Breal,Rimag=Bima,Preal=Rreal,Pimag=Pimag,RSS1=…RSSi=…RSS2 m =1.0,i∈[1,2 m ], where AXR and AXI are the real and imaginary parts of the intermediate vector (both initialized to zero vectors), Preal and Pimag are the real and imaginary parts of the search direction vector, Rreal and Rimag are the real and imaginary parts of the residual vector, RSSi is the error corresponding to the i-th pattern matrix, and the solution vectors XR and XI are the real and imaginary parts of the solution vector (both initialized to zero vectors).

[0119] After initialization, we can further search for the intermediate vectors AXR and AXI, the direction vectors Preal and Pimag, the residual vectors Rreal and Rimag, the residual norm VNRM, the error RSS1,…,RSSi…,RSS2 m , use the SIMD extension component to iteratively update the solution vectors XR and XI, and finally output the iterative solution vectors XR and XI. Then decompose the solution vectors XR and XI (decompose into corresponding 2 m The solutions of the pattern matrices, and correspondingly merge the real and imaginary parts of the solutions of each pattern matrix), we can get 2 m The solution vector of the pattern matrix.

[0120] Exemplarily, the parallel iterative solution process (corresponding to Figure 6 The process of loop iteration is as follows:

[0121] S1: Determine whether the iteration termination condition is currently met.

[0122] The iteration termination condition here is designed as follows: the number of iterations iter reaches the maximum number of iterations IMAX, or the error RSS of each pattern matrix at the current number of iterations is not greater than the convergence tolerance TOL.

[0123] S2: If the iteration termination condition is met, the iteration ends and the iterative solution vectors XR and XI are output.

[0124] S3: If the iteration termination condition is not met, start this iteration:

[0125] S31: Call the improved cnorm function to record the residual vector inner product value RNORM at the current number of iterations.

[0126] S32: Call the improved productAll function to calculate the intermediate vectors AXR and AXI at the current number of iterations.

[0127] S33: Calculate the update step size α at the current number of iterations based on the residual vector inner product value RNORM and the intermediate vectors AXR and AXI at the current number of iterations.

[0128] S34: Based on the update step size α and the search direction vectors Preal and Pimag at the current iteration number, the solution vectors XR and XI at the current iteration number are updated.

[0129] S35: Based on the update step size α and the intermediate vectors AXR and AXI at the current number of iterations, the residual vectors Rreal and Rimag are updated to obtain the residual vectors Rreal′ and Rimag′ at the current number of iterations.

[0130] S36: Calculate the conjugate parameter γ based on the residual vector inner product value RNORM at the current iteration number, and the residual vectors Rreal′ and Rimag′ at the current iteration number.

[0131] S37: Update the search direction vectors Preal and Pimag based on the conjugate parameter γ.

[0132] S38: The number of iterations is increased by 1, and the residual norm VNRN2 after the iteration is calculated. Based on the initialized residual norm VNRN and the residual norm VNRM2 after the iteration, the error RSS1,…,RSSi…,RSS2 under the current number of iterations is calculated. m , and then return to S1.

[0133] In order to introduce this solution in detail, this iteration process (i.e., the iteration process under the current iteration number) is described in detail here.

[0134] First, in S31, the improved cnorm function needs to be called to record the residual vector inner product value RNORM at the current number of iterations.

[0135] Initialize vectors sum1 and sum2.

[0136] SIMD extension unit reads residual vector Rreal2 m data, read 2 of the residual vector Rimag m data, where the SIMD extension unit reads data in steps of 2 m , Rreal is the real part of the residual vector R before updating at the current number of iterations, and Rimag is the imaginary part of the residual vector R before updating at the current number of iterations.

[0137] Then, calculate 2 according to the following formulam The residual vector inner product value RNORM:

[0138] RNORM=R T R, (1)

[0139] According to the improved cnorm function, we have:

[0140] RNORMreal=Rreal T Rreal-Rimag T Rimag, (2)

[0141] RNORMimag=Rreal T Rimag+Rimag T Real, (3)

[0142] Among them, RNORMreal is the real part of the residual vector inner product value RNORM, RNORMimag is the imaginary part of the residual vector inner product value RNORM, and Rreal T is the transpose of Rreal, Rimag T is the transpose of Rimag.

[0143] Then, CCR[0]=sum1[0], CCR[1]=sum1[1], ..., CCR[k]=sum1[k], ..., CCR[2 m -1]=sum1[2 m -1], save 2 m The real part of the inner product value of the residual vector RNORMreal is obtained by CCR[[2 m ]=sum2[0]、CCR[[2 m +1] = sum2[1], ..., CCR[2 m + k] = sum2[k], ..., CCR[2 m+1 -1]=sum2[2 m -1] Save 2 m The imaginary part of the inner product value of the residual vector RNORMimag.

[0144] Due to modern processor architectures (mainstream ones include x86 architecture and arm architecture, ARM NEON is the SIMD instruction set under the ARM platform, and the x86 platform has MMX, SSE, AVX and other instruction sets, and some other architectures also have some SIMD instruction sets), it can support the execution of multiple instructions simultaneously through its SIMD extension components (vectorization components). In this embodiment, NEON is used as an example to introduce the logic and pseudocode of the entire algorithm, which should not be regarded as a limitation of this application.

[0145] Taking the case of two pattern matrices as an example, the improved cnorm function (that is, the inner product function of the improved vector) is as follows:

[0146]

[0147]

[0148] As shown above, the vdupq_n_f64, vld1q_f64, and vmlaq_f64 functions of the NEON register are used for parallel calculation to realize the dot product calculation of two sets of vectors at a time. In this way, the calculation of the residual vector inner product value RNORM under the current number of iterations can be completed.

[0149] Secondly, in S32, the improved productAll function is called to calculate the intermediate vectors AXR and AXI at the current number of iterations.

[0150] The SIMD extension reads the search direction vectors Preal and Pimag, and fuses the new_CSR format quadruple: new_ValuesR array, new_ValuesI array, rowPtr array, and colIndices array. The SIMD extension reads data in steps of 2. m .

[0151] Then, the intermediate vectors AXR and AXI at the current iteration number are calculated according to the following formula:

[0152] AX=A*P, (4)

[0153] According to the improved complex multiplication function, we have:

[0154] AXR=AR*Preal-AI*Pimag, (5)

[0155] AXI=AR*Pimag+AI*Preal, (6)

[0156] Among them, AXR is the real part of the intermediate vector AX, AR is the real part of the pattern matrix, which is determined by reading the new_ValuesR array, rowPtr array and colIndices array through the SIMD extension component, Preal is the real part of the search direction vector P, AXI is the imaginary part of the intermediate vector AI, AI is the imaginary part of the pattern matrix, which is determined by reading the new_ValuesI array, rowPtr array and colIndices array through the SIMD extension component, and Pimag is the imaginary part of the search direction vector P.

[0157] Taking the case of NEON iteratively solving two pattern matrices as an example, the improved productAll function (i.e. the improved matrix multiplication vector function) is as follows:

[0158]

[0159]

[0160] When implementing the steps of matrix multiplication by vector, only selecting the non-zero value position for multiplication can greatly save time. For example, the first for loop here is used to control the i-th row of the matrix, that is, to select the i-th row to multiply the vector; the second for loop is used to control the non-zero value at the rowsPtr[i] position in the i-th row and the value at the colsIdx[j] position of the vector to perform multiplication. In this way, the calculation of the intermediate vectors AXR and AXI at the current number of iterations can be completed.

[0161] In S33, it is necessary to calculate the update step size α at the current number of iterations based on the residual vector inner product value RNORM at the current number of iterations and the intermediate vectors AXR and AXI.

[0162] The SIMD extension component reads the residual vector inner product values ​​RNORMreal and RNORMimag, the search direction vectors Preal and Pimag, and the intermediate vectors AXR and AXI at the current iteration number. The SIMD extension component reads data with a step size of 2. m .

[0163] Then, the update step size α at the current iteration number can be calculated according to the following formula:

[0164]

[0165] Further:

[0166]

[0167] Among them, αr is the real part of the update step size α, and αi is the imaginary part of the update step size α.

[0168] Here we need to use improved complex multiplication function and improved complex division function to implement it.

[0169] Taking the case of NEON iteratively solving two pattern matrices as an example, the improved complex multiplication function (ComMul_NEON function) is as follows:

[0170]

[0171] In the multiplication of complex numbers, NEON is used to calculate the multiplication of two sets of complex numbers at the same time to achieve parallel calculation. The input is also virtual, because the multiplication of complex numbers will be used in more than one place in the entire iteration process. Each time the corresponding input is made, the parallel multiplication calculation of complex numbers can be realized. The complex multiplication formula can be summarized as follows: for complex numbers Z1 = (a + bi), Z2 = (c + di), we have:

[0172] Z1*Z2=(a+bi)(c+di)=(ac-bd)+(ad+cb)i, (10)

[0173] Taking the case of NEON iteratively solving two pattern matrices as an example, the improved complex division function (ComDivi_NEON function) is as follows:

[0174]

[0175]

[0176] In the complex number division, NEON is used to calculate the division of two sets of complex numbers at the same time to achieve parallel calculation. The input is also virtual, because the entire iteration process will use complex number division in more than one place, and each time the corresponding input is made, the parallel division calculation of the complex number can be realized. The complex number division formula can be summarized as follows: for complex numbers Z1 = (a + bi), Z2 = (c + di), we have:

[0177]

[0178] Thus, the calculation of the update step length α under the current number of iterations can be completed. Note that the update step length α here is equal to the step length of the SIMD extension component reading data, which is 2. m They are two different concepts and should not be confused.

[0179] In S34, the solution vectors XR and XI at the current iteration number can be updated based on the update step size α and the search direction vectors Preal and Pimag at the current iteration number. Note that the search direction vectors Preal and Pimag here are not updated at the current iteration number, and the data updated in the last iteration are used, and the update order of the search direction vectors Preal and Pimag is later.

[0180] The SIMD extension component reads the update step size α and the search direction vectors Preal and Pimag at the current iteration number, where the step size of the SIMD extension component reading data is 2 m .

[0181] Then, the solution vector at the current iteration number is calculated and updated according to the following formula:

[0182] X=X+α*P, (12)

[0183] According to the improved complex multiplication function, we have:

[0184] XR′=XR+αr*Preal-αi*Preal, (13)

[0185] XI′=XI+αr*Pimag+αi*Preal, (14)

[0186] Among them, XR′ is the real part of the solution vector X after update under the current iteration number, XR is the real part of the solution vector X before update under the current iteration number, that is, the real part of the solution vector X updated in the last iteration, αr is the real part of the update step α, Preal is the real part of the search direction vector P, XI′ is the imaginary part of the solution vector X after update under the current iteration number, XI is the imaginary part of the solution vector X before update under the current iteration number, that is, the imaginary part of the solution vector X updated in the last iteration, αi is the imaginary part of the update step α, and Pimag is the imaginary part of the search direction vector P.

[0187] When updating the solution vector at the current number of iterations, the improved complex multiplication function (ComMul_NEON function) is used, which can be referred to in the previous article and will not be repeated here. In this way, the solution vector at the current number of iterations can be updated.

[0188] In S35, the residual vectors Rreal and Rimag can be updated based on the update step size α and the intermediate vectors AXR and AXI at the current iteration number to obtain the residual vectors Rreal′ and Rimag′ at the current iteration number. Note that the update step size α and the intermediate vectors AXR and AXI here are calculated at the current iteration number.

[0189] Based on this, the SIMD extension component reads the update step size α and the intermediate vectors AXR and AXI under the current iteration number, where the step size of the SIMD extension component reading data is 2 m .

[0190] Then, the residual vector under the current iteration number can be calculated and updated according to the following formula:

[0191] R′=R-α*AX, (15)

[0192] According to the improved complex multiplication function, we have:

[0193] Rreal′=Rreal-αr*AXR+αi*AXI, (16)

[0194] Rimag′=Rimag-αr*AXI-αi*AXR, (17)

[0195] Among them, Rreal′ is the real part of the residual vector R after update at the current iteration number, Rreal is the real part of the residual vector R before update at the current iteration number, that is, the real part of the residual vector R updated in the previous iteration, αr is the real part of the update step α, AXR is the real part of the intermediate vector AX at the current iteration number, Rimag′ is the imaginary part of the residual vector R after update at the current iteration number, Rimag is the imaginary part of the residual vector R before update at the current iteration number, that is, the imaginary part of the residual vector R updated in the previous iteration, αi is the imaginary part of the update step α, and AXI is the imaginary part of the intermediate vector AX at the current iteration number.

[0196] Similarly, when updating the residual vector at the current number of iterations, the improved complex multiplication function (ComMul_NEON function) is also used, which can be referred to in the previous article and will not be repeated here. In this way, the residual vector at the current number of iterations can be updated.

[0197] In S36, the conjugate parameter γ may be calculated based on the residual vector inner product value RNORM at the current iteration number, and the residual vectors Rreal′ and Rimag′ at the current iteration number.

[0198] The SIMD extension component reads the residual vector inner product value RNORM at the current iteration number, the residual vector Rreal′ and Rimag′ at the current iteration number, where the step length of the SIMD extension component reading data is 2 m .

[0199] Then, the conjugate parameter γ is calculated according to the following formula:

[0200]

[0201] Among them, γ is the conjugate parameter, and the conjugate parameter γ includes the real part γr and the imaginary part γi. Based on the improved complex multiplication function and the improved complex division function, we have:

[0202]

[0203] Among them, Rreal′ is the real part of the residual vector R′ after updating at the current number of iterations, Rreal′ T is the transpose of Rreal′, Rimag′ is the imaginary part of the residual vector R′ updated at the current number of iterations, and Rimag′ T is the transpose of Rimag′, Rreal is the real part of the residual vector R before updating at the current iteration number, that is, the real part of the residual vector R updated in the last iteration, Rreal Tis the transpose of Rreal, Rimag is the imaginary part of the residual vector R before update at the current iteration number, that is, the imaginary part of the residual vector R updated in the last iteration, Rimag T is the transpose of Rimag.

[0204] When calculating the conjugate parameter γ at the current number of iterations, the improved complex multiplication function (ComMul_NEON function) and the improved complex division function (ComDivi_NEON function) are also used, which can be referred to in the previous article and will not be repeated here. In this way, the conjugate parameter γ at the current number of iterations can be calculated.

[0205] In S37, the search direction vectors Preal and Pimag may be updated based on the conjugate parameter γ.

[0206] The SIMD extension component reads the search direction vectors Preal and Rimag, the real part γr and imaginary part γi of the conjugate parameter γ at the current iteration number, and the residual vectors Rreal′ and Rimag′ at the current iteration number (note that they are the updated residual vectors at the current iteration number). The step size of the SIMD extension component reading data is 2 m .

[0207] Then, the updated search direction vectors Preal and Pimag are calculated according to the following formula:

[0208] P′=R′+γP, (21)

[0209] According to the improved complex multiplication function, we have:

[0210] Rreal′=Rreal′+γr*Preal-γi*Pimag, (22)

[0211] Pimag′=Rimag′+γr*Pimag+γl*Preal, (23)

[0212] Among them, Preal′ is the real part of the search direction vector P′ after update at the current iteration number, Preal is the real part of the search direction vector P before update at the current iteration number, Pimag′ is the imaginary part of the search direction vector P′ after update at the current iteration number, and Pimag is the imaginary part of the search direction vector P before update at the current iteration number.

[0213] When updating the search direction vector at the current iteration number, the improved complex multiplication function (ComMul_NEON function) is also used, which can be referred to in the previous article and will not be repeated here. In this way, the search direction vector at the current iteration number can be updated.

[0214] In S38, after the above iterative update is completed, the number of iterations is increased by 1, that is, iter′=iter+1. Then, the residual norm VNRM2 after the iteration can be calculated.

[0215] The SIMD extension component reads the residual vectors Rreal′ and Rimag′ at the current number of iterations, where the step length of the SIMD extension component reading data is 2 m .

[0216] The residual norm VNRM2 of each mode matrix after this iteration is calculated according to the following formula:

[0217]

[0218] Where VNRM2 is the residual norm calculated by the pattern matrix at the current number of iterations, R′ is the updated residual vector at the current number of iterations, n is the dimension of the number of rows or columns in the pattern matrix, and Rreal′ j is the jth value of Rreal′, Rimag′ j is the jth value of Rimag′.

[0219] Taking the case of NEON iteratively solving two pattern matrices as an example, the improved vnorm function is used to calculate the residual norm VNRM2 of each pattern matrix after this iteration. The improved vnorm function (vnorm_CSR_NEON() function) is as follows:

[0220]

[0221] As mentioned above, the vdupq_n_f64() function of the NEON register is used to initialize the vector, and the vld1q_f64() function obtains two values ​​of the array to construct the vector and load it into Z1PF. Z1PF is a variable of type Float64x2_t, which is used to store two 64-bit float type variables. The vmlaq_f64() function is used to perform multiplication and addition operations to achieve parallel calculation of two matrices. Thus, the residual norm VNRM2 after iteration corresponding to each pattern matrix can be calculated.

[0222] After that, the errors RSS1,…,RSSi…,RSS2 under the current number of iterations can be further calculated m .

[0223] The SIMD extension component reads the initialized residual norm VNRM and the iterated residual norm VNRM2, where the step length of the SIMD extension component reading data is 2 m .

[0224] Then, the error of each pattern matrix at the current iteration number is calculated according to the following formula:

[0225]

[0226] Among them, RSSi is the error of the i-th pattern matrix at the current number of iterations, i∈[1,2 m ].

[0227] The above is all the processes of an iteration cycle. After completing the iteration update of this cycle, we can return to S1 to determine whether the iteration termination condition is met. The iteration termination condition is designed as follows: the number of iterations iter reaches the maximum number of iterations IMAX, or the error RSS of each pattern matrix at the current number of iterations is not greater than the convergence tolerance TOL.

[0228] In order to verify the effect of this solution, a test case is provided here (this embodiment takes NEON register as an example, as a SIMD extension component, and runs the iterative process in the data processing method for SIMD parallel iterative solution):

[0229] The test case is two complex sparse matrices of size 122036*122036, where the non-zero value of each matrix is ​​8168600. The comparison method is to perform CG iterative solution on the two matrices without using NEON registers, and record the solution time. The following is the comparison result of the two methods:

[0230] Unit\Method NEON parallel iterative solution Iterative solution NEON efficiency improvement Time(s) 2696.53 8087.68 66.66％

[0231] It can be seen that the solution time of this scheme and the comparison scheme is 2696.53 seconds and 8087.68 seconds respectively. The parallel iterative solution efficiency of this scheme is improved by 66.66%. If more matrices are solved, the acceleration performance will be more obvious.

[0232] In summary, the present application embodiment provides a data processing method for SIMD parallel iterative solution, using SIMD extension components, design 2 m A parallel iterative solution for the pattern matrix is ​​proposed to achieve accelerated solution. Specifically, the data to be processed (including 2 m pattern matrix, m∈Z +, each pattern matrix has the same non-zero element distribution pattern and is a complex matrix); store each pattern matrix in the real number new_CSR format, where the real number new_CSR format contains the ValuesR array, the ValuesI array, the rowPtr array, and the colIndices array. The ValuesR array is used to store the real part of the non-zero elements of the pattern matrix, the ValuesI array is used to store the imaginary part of the non-zero elements of the pattern matrix, the rowPtr array is used to store the position of the first non-zero element of each row in the pattern matrix in the values ​​array, and the colIndices array is used to store the column index of each non-zero element in the pattern matrix; for 2 m The real number new_CSR format of the pattern matrix is ​​vectorized to obtain the vectorized fused new_CSR format, where the fused new_CSR format includes the new_ValuesR array, the new_ValuesI array, the rowPtr array and the colIndices array, and the new_ValuesR array includes several real part arrays, each of which is used to store 2 m The real part of the non-zero elements at the same position of the pattern matrix, the new_ValuesI array contains the corresponding number of imaginary arrays, each imaginary array is used to store 2 m The imaginary part of the non-zero elements of the same position of the pattern matrix; use the SIMD extension component to read the vectorized fusion new_CSR format data, perform parallel iterative solution, and output the iterative solution result. This solution improves the storage format of CSR and vectorizes the corresponding data of multiple pattern matrices so that the SIMD extension component (such as NEON register) can read data and execute SIMD instructions (such as NEON instructions) in parallel, improves the relevant iterative solution algorithm, and thus achieves 2 m The parallel iterative solution of the pattern matrix greatly improves the iterative solution efficiency.

[0233] By storing each pattern matrix in a complex CSR format, the complex CSR format contains the Values ​​array, rowPtr array, and colIndices array. The Values ​​array is used to store the non-zero elements of the pattern matrix, the rowPtr array is used to store the position of the first non-zero element of each row in the pattern matrix in the values ​​array, and the colIndices array is used to store the column index of each non-zero element in the pattern matrix; for the complex CSR format of each pattern matrix: split the complex non-zero elements of the Values ​​array in the complex CSR format into the real part of the non-zero elements and the imaginary part of the non-zero elements, store the real part of the non-zero elements in sequence to obtain the ValuesR array, store the imaginary part of the non-zero elements in sequence to obtain the ValuesI array, and the rowPtr array and colIndices array remain unchanged. In this way, the existing complex CSR format can be used for further splitting and processing to obtain a new real number new_CSR format, which can make full use of existing technologies and be applicable to a wide range of application scenarios. And 2 m The real parts of the non-zero elements of the ValuesR array at the same position in the new_CSR format of the pattern matrix are combined into a 2-bit array in the order of the pattern matrix. m The real part array of the real part of the non-zero elements is obtained, and the new_ValuesR array containing several real part arrays is obtained; 2 m The imaginary parts of the non-zero elements in the ValuesI array at the same position in the real number new_CSR format of the pattern matrix are combined into a 2 m The imaginary part array of the imaginary part of the non-zero elements is obtained to obtain the new_ValuesR array containing several imaginary part arrays; the rowPtr array remains unchanged; the colIndices array remains unchanged. Since the use of vector registers to accelerate the program requires the data in the memory to be read into it in advance, vectorization is performed based on the real number new_CSR format to obtain the fused new_CSR format, which can achieve vectorization and make the fused new_CSR format adapt to the SIMD extension component. Data alignment can significantly improve the performance of the SIMD extension component, realize parallel iterative solution calculation, effectively improve the calculation speed of parallel processing, greatly improve calculation efficiency, and shorten iterative solution time.

[0234] The parallel iterative solution provided in this embodiment improves the iterative solution algorithm so that the entire iterative solution algorithm can be applied to 2 mThe parallel iterative solution calculation of the pattern matrix is ​​designed, and the corresponding improved formula is designed to realize parallel iterative update. At the same time, the operation step size of the SIMD extension component (such as reading data) is adjusted to ensure the stable operation of the parallel iterative solution. The test case uses this solution (taking the parallel accelerated iterative solution of NEON register as an example) to iteratively solve two complex sparse matrices of size 122036*122036, where the non-zero value of each matrix is ​​8168600. The comparison method is not to use the SIMD extension component to accelerate the iterative solution, and perform CG iterative solution on the two matrices respectively. The solution time is 2696.53 seconds and 8087.68 seconds respectively, and the iterative solution efficiency is improved by 66.66%. If more matrices are solved, the acceleration performance will be more obvious.

[0235] The above description is only an embodiment of the present application and is not intended to limit the protection scope of the present application. For those skilled in the art, the present application may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A data processing method for SIMD parallel iterative solution, characterized in that: include: Get the data to be processed, where the data to be processed includes 2 m pattern matrix, m∈Z + , each pattern matrix has the same non-zero element distribution pattern and is a complex matrix; Store each pattern matrix in the real number new_CSR format, where the real number new_CSR format contains the ValuesR array, the ValuesI array, the rowPtr array, and the colIndices array. The ValuesR array is used to store the real part of the non-zero elements of the pattern matrix, the ValuesI array is used to store the imaginary part of the non-zero elements of the pattern matrix, the rowPtr array is used to store the position of the first non-zero element of each row in the pattern matrix in the values ​​array, and the colIndices array is used to store the column index of each non-zero element in the pattern matrix. 2 m The real number new_CSR format of the pattern matrix is ​​vectorized to obtain the vectorized fused new_CSR format, where the fused new_CSR format includes the new_ValuesR array, the new_ValuesI array, the rowPtr array and the colIndices array, and the new_ValuesR array includes several real part arrays, each of which is used to store 2 m The real part of the non-zero elements at the same position of the pattern matrix, the new_ValuesI array contains the corresponding number of imaginary arrays, each imaginary array is used to store 2 m The imaginary parts of the non-zero elements of the pattern matrices at the same position; The SIMD extension component is used to read the vectorized fused new_CSR format data, perform parallel iterative solutions, and output the iterative solution results.

2. The data processing method for SIMD parallel iterative solution according to claim 1, characterized in that: Store each pattern matrix in real new_CSR format, consisting of: Store each pattern matrix in a complex CSR format, where the complex CSR format contains a Values ​​array, a rowPtr array, and a colIndices array. The Values ​​array is used to store the non-zero elements of the pattern matrix. The rowPtr array is used to store the position of the first non-zero element of each row in the pattern matrix in the values ​​array. The colIndices array is used to store the column index of each non-zero element in the pattern matrix. For the complex CSR format of each pattern matrix: split the non-zero elements in the complex form of the Values ​​array in the complex CSR format into the real part of the non-zero elements and the imaginary part of the non-zero elements, store the real part of the non-zero elements in sequence to obtain the ValuesR array, store the imaginary part of the non-zero elements in sequence to obtain the ValuesI array, and the rowPtr array and colIndices array remain unchanged.

3. The data processing method for SIMD parallel iterative solution according to claim 1, characterized in that: 2 m The real number new_CSR format of the pattern matrix is ​​vectorized to obtain the vectorized fusion new_CSR format, including: 2 m The real parts of the non-zero elements of the ValuesR array at the same position in the new_CSR format of the pattern matrix are combined into a 2-bit array in the order of the pattern matrix. m The real part array of the real part of non-zero elements is obtained, and the new_ValuesR array containing several real part arrays is obtained; 2 m The imaginary parts of the non-zero elements in the ValuesI array at the same position in the real number new_CSR format of the pattern matrix are combined into a 2 m The imaginary part array of the non-zero elements of the imaginary part is obtained, and the new_ValuesR array containing several imaginary part arrays is obtained; The rowPtr array remains unchanged; The colIndices array remains unchanged.

4. The data processing method for SIMD parallel iterative solution according to claim 1, characterized in that: Use SIMD extension components to read vectorized fused new_CSR format data and perform parallel iterative solutions, including: Using SIMD extensions to read 2 m The fused new_CSR format quaternion of the pattern matrix, the fused new_CSR format right value Breal and Bimag, the dimension n of the pattern matrix, the maximum number of iterations IMAX, and the convergence tolerance TOL, where the fused new_CSR format quaternion is 2 m The matrix A in the pattern matrix Ax=b is stored in the real number new_CSR format and then vectorized, including the new_ValuesR array, new_ValuesI array, rowPtr array and colIndices array, and the right value of the fused new_CSR format is 2 m The right-hand side vector b in the pattern matrix Ax=b is obtained by storing it in the real number new_CSR format and then quantizing it, including Breal and Bimag, Breal is the real part of the right value, and Bimag is the imaginary part of the right value; Initialize intermediate vectors AXR and AXI, search direction vectors Preal and Pimag, residual vectors Rreal and Rimag, solution vectors XR and XI, residual norm VNRM, errors RSS1,…,RSSi…,RSS2 m , and assign values: Rreal=Bral,Rimag=Bima,Preal=Rreal,Piamg=Rimag,RSS1=…RSSi=…RSS2 m =1.0,i∈[1,2 m ], where AXR and AXI are the real and imaginary parts of the intermediate vector, Preal and Pimag are the real and imaginary parts of the search direction vector, Rreal and Rimag are the real and imaginary parts of the residual vector, and RSSi is the error corresponding to the i-th pattern matrix; Based on the intermediate vectors AXR and AXI, the search direction vectors Preal and Pimag, the residual vectors Rreal and Rimag, the residual norm VNRM, the errors RSS1,…,RSSi…,RSS2 m , using the SIMD extension component to iteratively update the solution vectors XR and XI, and finally output the iteratively solved solution vectors XR and XI; Decompose the solution vectors XR and XI to obtain 2 m The solution vector of the pattern matrix.

5. The data processing method for SIMD parallel iterative solution according to claim 4, characterized in that: Based on the intermediate vectors AXR and AXI, the search direction vectors Preal and Pimag, the residual vectors Rreal and Rimag, the residual norm VNRM, the errors RSS1,…,RSSi…,RSS2 m , using the SIMD extension component to iteratively update the solution vectors XR and XI, and finally output the iteratively solved solution vectors XR and XI, including: S1: Determine whether the iteration termination condition is met; S2: If the iteration termination condition is met, the iteration ends and the iterative solution vectors XR and XI are output; S3: If the iteration termination condition is not met, start this iteration: S31: Call the improved cnorm function to record the residual vector inner product value RNORM at the current number of iterations; S32: Call the improved productAll function to calculate the intermediate vectors AXR and AXI at the current number of iterations; S33: Calculate the update step size α at the current number of iterations based on the residual vector inner product value RNORM and the intermediate vectors AXR and AXI at the current number of iterations; S34: based on the update step size α and the search direction vectors Preal and Pimag at the current iteration number, update the solution vectors XR and XI at the current iteration number; S35: based on the update step size α and the intermediate vectors AXR and AXI at the current iteration number, the residual vectors Rreal and Rimag are updated to obtain the residual vectors Rreal′ and Rimag′ at the current iteration number; S36: Calculate the conjugate parameter γ based on the residual vector inner product value RNORM at the current iteration number, the residual vectors Rreal′ and Rimag′ at the current iteration number; S37: Update the search direction vectors Preal and Pimag based on the conjugate parameter γ; S38: The number of iterations is increased by 1, and the residual norm VNRN2 after the iteration is calculated. Based on the initial residual norm VNRM and the residual norm VNRM2 after the iteration, the error RSS1,…,RSSi…,RSS2 under the current number of iterations is calculated. m , and then return to S1.

6. The data processing method for SIMD parallel iterative solution according to claim 5, characterized in that: Call the improved cnorm function to record the residual vector inner product value RNORM at the current number of iterations, including: Initialize vectors sum1 and sum2; SIMD extension unit reads residual vector Rreal2 m data, read 2 of the residual vector Rimag m data, where the SIMD extension unit reads data in steps of 2 m , Rreal is the real part of the residual vector R before updating at the current number of iterations, and Rimag is the imaginary part of the residual vector R before updating at the current number of iterations; Calculate according to the following formula2 m The residual vector inner product value RNORM: RNORM=RTR, RNORMreal=RrealTRreal-RimagTRimag, RNORMimag=RrealTRimag+RimagTRreal, Among them, RNORMreal is the real part of the residual vector inner product value RNORM, RNORMimag is the imaginary part of the residual vector inner product value RNORM, RrealT is the transpose of Rreal, and RimagT is the transpose of Rimag; By CCR[0]=sum1[0], CCR[1]=sum1[1], ..., CCR[k]=sum1[k], ..., CCR[2 m -1]=sum1[2 m -1] Save 2 m The real part of the inner product value of the residual vector RNORMreal is obtained by CCR[[2 m ]=sum2[0]、CCR[[2 m +1] = sum2[1], ..., CCR[2 m + k] = sum2[k], ..., CCR[2 m+1 -1]=sum2[2 m -1] Save 2 m The imaginary part of the inner product value of the residual vector RNORMimag.

7. The data processing method for SIMD parallel iterative solution according to claim 5, characterized in that: Call the improved productAll function to calculate the intermediate vectors AXR and AXI at the current number of iterations, including: The SIMD extension reads the search direction vectors Preal and Pimag, and fuses the new_CSR format quadruple: new_ValuesR array, new_ValuesI array, rowPtr array, and colIndices array. The SIMD extension reads data in steps of 2. m ; The intermediate vectors AXR and AXI at the current iteration number are calculated according to the following formula: AXR=AR*Preal-AI*Pimag, AXI=AR*Pimag+AI*Preal, Among them, AXR is the real part of the intermediate vector AX, AR is the real part of the pattern matrix, which is determined by reading the new_ValuesR array, rowPtr array and colIndices array through the SIMD extension component, Preal is the real part of the search direction vector P, AXI is the imaginary part of the intermediate vector AI, AI is the imaginary part of the pattern matrix, which is determined by reading the new_ValuesI array, rowPtr array and colIndices array through the SIMD extension component, and Pimag is the imaginary part of the search direction vector P.

8. The data processing method for SIMD parallel iterative solution according to claim 6, characterized in that: Based on the residual vector inner product value RNORM and the intermediate vectors AXR and AXI at the current number of iterations, the update step size α at the current number of iterations is calculated, including: The SIMD extension component reads the residual vector inner product values ​​RNORMreal and RNORMimag, the search direction vectors Preal and Pimag, and the intermediate vectors AXR and AXI at the current iteration number. The SIMD extension component reads data with a step size of 2. m ; The update step size α at the current number of iterations is calculated according to the following formula: Further: Among them, αr is the real part of the update step size α, and αi is the imaginary part of the update step size α.

9. The data processing method for SIMD parallel iterative solution according to claim 8, characterized in that: Based on the update step size α and the search direction vectors Preal and Pimag at the current iteration number, the solution vectors XR and XI at the current iteration number are updated, including: The SIMD extension component reads the update step size α and the search direction vectors Preal and Pimag at the current iteration number, where the step size of the SIMD extension component reading data is 2 m ; The solution vector for the current iteration number is calculated and updated according to the following formula: XR′=XR+αr+Preal-αi*Pimag, XI′=XI+αr*Pimag+αi*Preal, Among them, XR′ is the real part of the solution vector X after update under the current iteration number, XR is the real part of the solution vector X before update under the current iteration number, that is, the real part of the solution vector X updated in the last iteration, αr is the real part of the update step α, Preal is the real part of the search direction vector P, XI′ is the imaginary part of the solution vector X after update under the current iteration number, XI is the imaginary part of the solution vector X before update under the current iteration number, that is, the imaginary part of the solution vector X updated in the last iteration, αi is the imaginary part of the update step α, and Pimag is the imaginary part of the search direction vector P.

10. The data processing method for SIMD parallel iterative solution according to claim 8, characterized in that: Based on the update step size α and the intermediate vectors AXR and AXI at the current number of iterations, the residual vectors Rreal and Rimag are updated to obtain the residual vectors Rreal′ and Rimag′ at the current number of iterations, including: The SIMD extension unit reads the update step size α and the intermediate vectors AXR and AXI at the current iteration number, where the step size of the SIMD extension unit reading data is 2 m ; The residual vector for the current iteration number is calculated and updated according to the following formula: Rreal′=Rreal-αr*AXR+αi*AXI, Rimag′=Rimag-αr*AXI-αi*AXR, Among them, Rreal′ is the real part of the residual vector R after update at the current iteration number, Rreal is the real part of the residual vector R before update at the current iteration number, that is, the real part of the residual vector R updated in the previous iteration, αr is the real part of the update step α, AXR is the real part of the intermediate vector AX at the current iteration number, Rimag′ is the imaginary part of the residual vector R after update at the current iteration number, Rimag is the imaginary part of the residual vector R before update at the current iteration number, that is, the imaginary part of the residual vector R updated in the previous iteration, αi is the imaginary part of the update step α, and AXI is the imaginary part of the intermediate vector AX at the current iteration number.

Images