Data processing method for SIMD-oriented parallel iterative solution

Abstract

The present application provides a data processing method for parallel iterative solution oriented to SIMD. By obtaining data to be processed (including 2<supgt;m< / supgt; complex form pattern matrices), each pattern matrix is stored in a real number new_CSR format (including a ValuesR array, a ValuesI array, a rowPtr array, and a colIndices array), and vectorization is performed to obtain a vectorized fused new_CSR format (including a new_ValuesR array, a new_ValuesI array, a rowPtr array, and a colIndices array). The SIMD extension component is used to read the vectorized fused new_CSR format data for parallel iterative solution, and the iterative solution result is output. In this way, parallel iterative solution calculation can be performed through the SIMD extension component in the SIMD architecture to achieve accelerated solution.

Description

Technical Field

[0001] This application relates to the technical field of data processing, and more specifically, to a data processing method for parallel iterative solution oriented to SIMD. Background Art

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

[0003] The iterative methods for linear equations are a class of methods for solving large sparse linear equations, which solve the equations by gradually approaching 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 introductions:

[0004] 1. Jacobi Iteration

[0005] Basic idea: Decompose the linear equation 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, approximate the solution step by step through the iterative formula x (k+1) = D -1 (b - (L + U)x (k) ).

[0006] Characteristics: Simple and easy to implement, but with a slow convergence rate, suitable for diagonally dominant matrices.

[0007] 2. Gauss - Seidel Iteration

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

[0009] Characteristics: Usually converges faster than the Jacobi iteration method, but with a slightly higher computational complexity.

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

[0011] Basic idea: Introduce a relaxation parameter ω on the basis of the Gauss - Seidel iteration method and accelerate convergence by adjusting ω. The iterative formula is:

[0012]

[0013] Feature: By selecting an appropriate ω, the convergence speed can be significantly accelerated. However, inappropriate selection may lead to divergence.

[0014] 4. Conjugate Gradient Method (CG)

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

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

[0017] 5. Biconjugate Gradient Method (BCG)

[0018] Basic idea: Applicable to non-symmetric matrices, dealing with non-symmetric cases by constructing two mutually conjugate sequences.

[0019] Feature: Applicable to non-symmetric matrices, but may encounter convergence problems, especially in the case of ill-conditioned matrices or matrices with complex eigenvalues.

[0020] 6. Biconjugate Gradient Stabilized Method (BiCGSTAB)

[0021] Basic idea: Based on the BCG method, an additional stabilization step is introduced to improve convergence and stability.

[0022] Feature: More stable than BCG, faster convergence speed, applicable to non-symmetric matrices.

[0023] 7. Generalized Minimal Residual Method (GMRES)

[0024] Basic idea: Gradually approaching the solution by minimizing the residual norm, applicable to non-symmetric matrices. A Krylov subspace is reconstructed at each step.

[0025] Feature: Good convergence, but relatively high computational complexity, and multiple vectors need to be stored.

[0026] 8. Preconditioned Conjugate Gradient Method (PCG)

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

[0028] Features: By selecting an appropriate preconditioner, the convergence rate can be significantly improved, which is applicable to symmetric positive definite matrices.

[0029] These methods each have their own advantages and disadvantages. The selection of an appropriate method depends on the specific nature of the problem, such as the symmetry, sparsity, and ill-conditioning degree of the matrix.

[0030] Modern processor architectures (the mainstream ones include the x86 architecture and the ARM architecture, etc. ARM NEON is a SIMD instruction set under the ARM platform, and the x86 platform has instruction sets such as MMX, SSE, and AVX. Some other architectures also have some SIMD instruction sets) can support the execution of multiple instructions simultaneously through their SIMD extension components (vectorization components). SIMD, that is, single instruction multiple data, means that a single operation instruction can execute multiple data streams, thereby improving the operation speed. The SIMD extension component refers to the dedicated hardware unit in the processor for executing 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 features.

[0031] 1. Intel's SIMD extension: Name: SSE (Streaming SIMD Extensions); Versions: SSE, SSE2, SSE3, SSSE3, SSE4.1, SSE4.2; Register width: 128 bits; Data types: single-precision floating-point numbers (32 bits), double-precision floating-point numbers (64 bits), integers (8 bits, 16 bits, 32 bits, 64 bits); Applications: widely used in floating-point operations, multimedia processing, and scientific computing; Registers: XMM0 to XMM15 (there may be 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 numbers (32 bits); Applications: mainly used in floating-point operations, especially in games and multimedia applications; Registers: MM0 to MM7.

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

[0034] Taking NEON as an example, the NEON technology is an advanced SIMD (Single Instruction Multiple Data) architecture for Arm Cortex-A series processors. It can accelerate multimedia and signal processing algorithms, such as video encoders / decoders, 2D / 3D graphics, games, audio and speech processing, image processing, telephony, and sound. NEON instructions perform "packed SIMD" processing, and registers are considered vectors of elements of the same data type. The supported data types are: signed / unsigned 8-bit, 16-bit, 32-bit, and 64-bit.

[0035] There are a total of 32 NEON vector registers in ARMv8, with a length of 128 bits. They can be used for processing scalar operands or vector operands. These registers can store multiple operands of the same data type, such as 2 double-type operands of 64 bits or 4 float-type operands of 32 bits. During the execution of NEON instructions, multiple operands on the same vector register are processed simultaneously, thus achieving parallel processing of different data. Accelerating a program using vector registers requires reading the data in memory into them in advance. Therefore, 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, the solution of systems of equations with the same non-zero element distribution positions is usually faced. For example, in the field of electromagnetic computing, when the model and mesh dissection remain unchanged, the non-zero element distribution positions of the systems of equations generated at different frequencies are the same. Since there are a large number of complex matrices (dense complex matrices, sparse complex matrices, etc.) that need to be linearly solved in the field of electromagnetic numerical analysis, but the current technology requires calling a solution function each time for solving, and there is no way to perform parallel solution of 2 or more linear equations.

[0037] Based on this, the inventors of the present application designed a parallel iterative solution scheme for SIMD, which can achieve parallel iterative solution of pattern matrices (large complex matrices with the same non-zero element distribution pattern) in various scenarios and 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 for SIMD, so as to perform parallel iterative solution calculation through SIMD extension components and achieve accelerated solution.

[0039] To achieve the above purpose, the embodiments of the present application are implemented as follows:

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

[0041] In combination with the first aspect, in the first possible implementation manner of the first aspect, each pattern matrix is stored in the real new_CSR format, including: storing each pattern matrix in the complex CSR format, where 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 in each row of the pattern matrix in the values array, and the colIndices array is used to store the column indices of each non-zero element in the pattern matrix. For the complex CSR format of each pattern matrix: split the complex non-zero elements in the Values array of the complex CSR format into the real part of the non-zero element and the imaginary part of the non-zero element, store the real part of the non-zero element in order to obtain the ValuesR array, store the imaginary part of the non-zero element in order to obtain the ValuesI array, and the rowPtr array and the colIndices array remain unchanged.

[0042] In combination with the first aspect, in the second possible implementation manner of the first aspect, vectorize the real new_CSR format of 2 m pattern matrices to obtain a vectorized fused new_CSR format, including: combining the real parts of the non-zero elements at the same position in the ValuesR arrays of the real new_CSR formats of 2 m pattern matrices in the order of the pattern matrices into a real part array containing 2 m real parts of non-zero elements to obtain a new_ValuesR array containing several real part arrays; combining the imaginary parts of the non-zero elements at the same position in the ValuesI arrays of the real new_CSR formats of 2 m pattern matrices in the order of the pattern matrices into an imaginary part array containing 2 m imaginary parts of non-zero elements to obtain a new_ValuesR array containing several imaginary part arrays; the rowPtr array remains unchanged; the colIndices array remains unchanged.

[0043] In combination with the first aspect, in the third possible implementation manner of the first aspect, use the SIMD extension component to read the vectorized fused new_CSR format data and perform parallel iterative solution, including: using the SIMD extension component to read the quadruple of the fused new_CSR format of 2 m pattern matrices, the right values Breal and Bimag of the fused new_CSR format, the dimension n of the pattern matrix, the maximum number of iterations IMAX, and the convergence tolerance TOL, where the quadruple of the fused new_CSR format is 2 mThe matrix A in the pattern matrix Ax = b is obtained by vectorizing after being stored in the real new_CSR format, and includes the new_ValuesR array, the new_ValuesI array, the rowPtr array, and the colIndices array. The right value of the fused new_CSR format is 2 m The right-hand vector b in the pattern matrix Ax = b is obtained by vectorizing after being stored in the real 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; initialize the intermediate vectors AXR and AXI, the search direction vectors Preal and Pimag, the residual vectors Rreal and Rimag, the solution vectors XR and XI, the residual norm VNRM, the errors RSS1,…, RSSi…, RSS2 m , and assign values: Rreal = Breal, 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 respectively, Preal and Pimag are the real and imaginary parts of the search direction vector respectively, Rreal and Rimag are the real and imaginary parts of the residual vector respectively, 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 iteratively solved solution vectors XR and XI; decompose the solution vectors XR and XI to obtain 2 m solution vectors of the pattern matrix

[0044] Combined with the third possible implementation manner of the first aspect, in the fourth possible implementation manner 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, the solution vectors XR and XI are iteratively updated using the SIMD extension component, and finally the iteratively solved solution vectors XR and XI are output, including: S1: Determine whether the current iteration termination condition is satisfied; S2: If the iteration termination condition is satisfied, end the iteration and output the iteratively solved solution vectors XR and XI; S3: If the iteration termination condition is not satisfied, start the current iteration: S31: Call the improved cnorm function to record the inner product value RNORM of the residual vector at the current iteration; S32: Call the improved productAll function to calculate the intermediate vectors AXR and AXI at the current iteration; S33: Calculate the update step α at the current iteration based on the inner product value RNORM of the residual vector and the intermediate vectors AXR and AXI at the current iteration; S34: Update the solution vectors XR and XI at the current iteration based on the update step α and the search direction vectors Preal and Pimag at the current iteration; S35: Update the residual vectors Rreal and Rimag based on the update step α and the intermediate vectors AXR and AXI at the current iteration to obtain the residual vectors Rreal' and Rimag' at the current iteration; S36: Calculate the conjugate parameter γ based on the inner product value RNORM of the residual vector and the residual vectors Rreal' and Rimag' at the current iteration; S37: Update the search direction vectors Preal and Pimag based on the conjugate parameter γ; S38: Increment the iteration count, calculate the residual norm VNRM2 after iteration, and calculate the error RSS1, …, RSSi …, RSS2 at the current iteration based on the initialized residual norm VNRM and the residual norm VNRM2 after iteration. m , and then return to S1.

[0045] Combined with the fourth possible implementation manner of the first aspect, in the fifth possible implementation manner of the first aspect, calling the improved cnorm function to record the inner product value RNORM of the residual vector at the current iteration includes: initializing the vectors sum1 and sum2; the SIMD extension component reads 2 m data from the residual vector Rreal and reads 2 m data from the residual vector Rimag, where the step size for the SIMD extension component to read data is 2 m , Rreal is the real part of the residual vector R before update at the current iteration, and Rimag is the imaginary part of the residual vector R before update at the current iteration; calculate 2 m inner product values RNORM of the residual vector according to the following formula:

[0046] RNORM = R T R,

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

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

[0049] Wherein, RNORMreal is the real part of the residual norm RNORM, RNORMimag is the imaginary part of the inner product value of the residual vector RNORM, Rreal T is the transpose of Rreal, Rimag T is the transpose of Rimag; through CCR[0] = sum1[0], CCR[1] = sum1[1], …, CCR[k] = sum1[k], …, CCR[2 m - 1] = sum1[2 m - 1] to save the real part RNORMreal of 2 m inner product values of the residual vector, and through 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] to save the imaginary part RNORMimag of 2 m inner product values of the residual vector.

[0050] Combined with the fourth possible implementation manner of the first aspect, in the sixth possible implementation manner of the first aspect, an improved productAll function is called to calculate the intermediate vectors AXR and AXI at the current iteration number, including: the SIMD extension component reads the search direction vectors Preal and Pimag, and fuses the quadruples in the new_CSR format: the new_ValuesR array, the new_ValuesI array, the rowPtr array, and the colIndices array, wherein the step size for the SIMD extension component to read 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 the SIMD extension component reading the new_ValuesR array, rowPtr array, and colIndices array. Preal is the real part of the search direction vector P. AXI is the imaginary part of the intermediate vector AI, and AI is the imaginary part of the pattern matrix, which is determined by the SIMD extension component reading the new_ValuesI array, rowPtr array, and colIndices array. Pimag is the imaginary part of the search direction vector P.

[0054] Combined with the fifth possible implementation manner of the first aspect, in the seventh possible implementation manner of the first aspect, based on the inner product value RNORM of the residual vector, and the intermediate vectors AXR and AXI at the current iteration number, calculate the update step size α at the current iteration number, including: the SIMD extension component reads the real part RNORMreal and imaginary part RNORMimag of the inner product value of the residual vector at the current iteration number, the real part Preal and imaginary part Pimag of the search direction vector, and the intermediate vectors AXR and AXI at the current iteration number. Among them, the step size for the SIMD extension component to read data is 2. m ;

[0055] Calculate the update step size α at the current iteration number according to the following formula:

[0056]

[0057] Furthermore:

[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] Combined with the seventh possible implementation manner of the first aspect, in the eighth possible implementation manner of the first aspect, based on the update step size α and the real part Preal and imaginary part Pimag of the search direction vector at the current iteration number, update the solution vectors XR and XI at the current iteration number, including: the SIMD extension component reads the update step size α and the real part Preal and imaginary part Pimag of the search direction vector at the current iteration number. Among them, the step size for the SIMD extension component to read data is 2. m ; Calculate the solution vectors for updating at the current iteration number according to the following formula:

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

[0062] XI′ = XI + ar * Pimag + αi * Preal,

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

[0064] Combined with the seventh possible implementation manner of the first aspect, in the ninth possible implementation manner of the first aspect, based on the update step size α, the intermediate vectors AXR and AXI at the current iteration, update the residual vectors Rreal and Rimag to obtain the residual vectors Rreal′ and Rimag′ at the current iteration, including: The SIMD extension component reads the update step size α, the intermediate vectors AXR and AXI at the current iteration, where the step size for the SIMD extension component to read data is 2 m ; Calculate and update the residual vector at the current iteration 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 updated residual vector R at the current iteration, Rreal is the real part of the residual vector R before update at the current iteration, that is, the real part of the residual vector R updated in the previous iteration, αr is the real part of the update step size α, AXR is the real part of the intermediate vector AX at the current iteration, Rimag′ is the imaginary part of the updated residual vector R at the current iteration, Rimag is the imaginary part of the residual vector R before update at the current iteration, that is, the imaginary part of the residual vector R updated in the previous iteration, α is the imaginary part of the update step size α, and AXI is the imaginary part of the intermediate vector AX at the current iteration.

[0068] Combined with the ninth possible implementation manner of the first aspect, in the tenth possible implementation manner of the first aspect, based on the inner product value RNORM of the residual vector at the current iteration, the residual vectors Rreal′ and Rimag′ at the current iteration, calculate the conjugate parameter γ, including: The SIMD extension component reads the inner product value RNORM of the residual vector at the current iteration, the residual vectors Rreal′ and Rimag′ at the current iteration, where the step size for the SIMD extension component to read data is 2 m ; Calculate the conjugate parameter γ according to the following formula:

[0069]

[0070] Among them, γ is the conjugate parameter. The conjugate parameter γ includes a real part γr and an imaginary part γi, and there is:

[0071]

[0072] Among them, Rreal′ is the real part of the updated residual vector R′ at the current iteration, Rreal′T is the transpose of Rreal′, Rimag′ is the imaginary part of the updated residual vector R′ at the current iteration, Rimag′ T is the transpose of Rimag′, Rreal is the real part of the residual vector R before update at the current iteration, that is, the real part of the residual vector R updated in the previous iteration, RrealT is the transpose of Rreal, Pimag is the imaginary part of the residual vector R before update at the current iteration, that is, the imaginary part of the residual vector R updated in the previous iteration, and RimagT is the transpose of Rimag.

[0073] Combined with the tenth possible implementation manner of the first aspect, in the eleventh possible implementation manner of the first aspect, updating the search direction vectors Preal and Pimag based on the conjugate parameter γ includes: 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 iteration, and the residual vectors Rreal′ and Rimag′ at the current iteration. Among them, the step size for the SIMD extension component to read 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 + γi * Preal,

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

[0077] Combined with the ninth possible implementation manner of the first aspect, in the twelfth possible implementation manner of the first aspect, calculate the residual norm VNRM2 after iteration, and calculate the error RSS1, …, RSSi, …, RSS2 at the current iteration number based on the initialized residual norm VNRM and the residual norm VNRM2 after iteration. m , including: the SIMD extension component reads the residual vectors Rreal′ and Rimag′ at the current iteration number, where the step size for the SIMD extension component to read data is 2. m ; calculate the residual norm VNRM2 of each mode matrix after this iteration according to the following formula:

[0078]

[0079] where VNRM2 is the residual norm calculated for the mode matrix at the current iteration number, R′ is the updated residual vector at the current iteration number, n is the dimension of the number of rows or the number of columns in the mode matrix, Rreal′ j is the j-th value of Rreal′, and Rimag′ j is the j-th value of Rimag′; the SIMD extension component reads the initialized residual norm VNRM and the residual norm VNRM2 after iteration, where the step size for the SIMD extension component to read data is 2. m ; calculate the error of each mode matrix at the current iteration number according to the following formula:

[0080]

[0081] where RSSi is the error of the i-th mode matrix at the current iteration number, and i ∈ [1, 2 m .

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

[0083] Beneficial effects:

[0084] 1. In the field of numerical analysis (such as the field of electromagnetic numerical analysis, the field of mechanical numerical analysis, etc.), there are a large number of complex sparse matrices that need to be linearly solved. However, each time a solution is required, a solution function needs to be called, and it is impossible to solve 2 or more linear equations in parallel. This solution uses the SIMD extension component to design a parallel iterative solution scheme for 2 m mode matrices to achieve accelerated solution. Specifically, by obtaining the data to be processed (including 2 mA 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 - valued new_CSR format, where the real - valued 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 parts of the non - zero elements of the pattern matrix, the ValuesI array is used to store the imaginary parts of the non - zero elements of the pattern matrix, the rowPtr array is used to store the position of the first non - zero element in each row of the pattern matrix in the values array, and the colIndices array is used to store the column indices of each non - zero element in the pattern matrix; for 2 m pattern matrices in the real - valued new_CSR format, vectorize them to obtain a vectorized fused new_CSR format. Among them, the fused new_CSR format includes the new_ValuesR array, the new_ValuesI array, the rowPtr array, and the colIndices array. The new_ValuesR array contains several real - part arrays, and each real - part array is used to store the real parts of the non - zero elements at the same position of 2 m pattern matrices. The new_ValuesI array contains the corresponding number of imaginary - part arrays, and each imaginary - part array is used to store the imaginary parts of the non - zero elements at the same position of 2 m pattern matrices; use the SIMD extension component to read the data in the vectorized fused new_CSR format and perform parallel iterative solution, and output the iterative solution result. This solution improves the storage format of CSR. By vectorizing and storing the corresponding data of multiple pattern matrices, the SIMD extension component (such as the NEON register) can parallelly read data and execute SIMD instructions (such as NEON instructions), and improve the relevant iterative solution algorithm, so as to achieve the parallel iterative solution of 2 m pattern matrices, greatly improving the iterative solution efficiency.

[0085] 2. This solution can store each pattern matrix in the complex CSR format. 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 in each row of the pattern matrix in the values array. The colIndices array is used to store the column indices of each non-zero element in the pattern matrix. For the complex CSR format of each pattern matrix: split the complex non-zero elements in the Values array of the complex CSR format into the real part of the non-zero element and the imaginary part of the non-zero element, store the real parts of the non-zero elements in order to obtain the ValuesR array, store the imaginary parts of the non-zero elements in order to obtain the ValuesI array, and keep the rowPtr array and the colIndices array unchanged. In this way, the existing complex CSR format can be used for further splitting processing to obtain a new real number new_CSR format, which can make full use of the existing technology and be applied to a wide range of application scenarios. And for the m real part of the non-zero elements at the same position in the ValuesR array of the real number new_CSR format of two m pattern matrices are combined in the order of the pattern matrices to form a real part array containing two m real parts of non-zero elements to obtain a new_ValuesR array containing several real part arrays; for the imaginary parts of the non-zero elements at the same position in the ValuesI array of the real number new_CSR format of two m pattern matrices are combined in the order of the pattern matrices to form an imaginary part array containing two

[0086] imaginary parts of non-zero elements to obtain a new_ValuesI array containing several imaginary part arrays; the rowPtr array remains unchanged; the colIndices array remains unchanged. Since accelerating the program using vector registers requires reading the data in memory into them in advance, therefore, vectorization is performed based on the real number new_CSR format to obtain a fused new_CSR format, which can achieve vectorization and make the fused new_CSR format adapt to the SIMD extension unit. The alignment of data can significantly improve the performance of the SIMD extension unit, realize parallel iterative solution calculation, effectively improve the calculation speed of parallel processing, greatly improve the calculation efficiency, and shorten the iterative solution time. m 3. The parallel iterative solution scheme of this solution improves the iterative solution algorithm so that the entire iterative solution algorithm can be applied to m 3. The parallel iterative solution scheme of this solution improves the iterative solution algorithm so that the entire iterative solution algorithm can be applied to mFor the calculation of parallel iterative solution of a pattern matrix, an improved formula is designed to achieve parallel iterative update. At the same time, the operation step (such as data reading) of the SIMD extension component is adjusted to ensure the stable operation of parallel iterative solution. Taking this solution (taking the parallel accelerated iterative solution using NEON registers as an example), test cases perform iterative solution on two complex sparse matrices of size 122036 * 122036, where the number of non-zero values in each matrix is 8168600. The comparison method is iterative solution without using SIMD extension component acceleration. The CG iterative solution is performed on the two matrices respectively, and the obtained solution times are 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 objects, features and advantages of the present application more obvious and understandable, the following specifically gives preferred embodiments and, in conjunction with the accompanying drawings, the detailed description is as follows. Brief Description of the Drawings

[0088] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the following will briefly introduce the drawings required to be used in the embodiments of the present application. It should be understood that the following drawings only show some embodiments of the present application, and therefore should not be regarded as limiting the scope. For those of ordinary skill in the art, other related drawings can be obtained based on these drawings without creative efforts.

[0089] Figure 1 It is a flowchart of a data processing method for parallel iterative solution oriented to SIMD provided by an embodiment of the present application.

[0090] Figure 2 It is a schematic diagram for storing a pattern matrix in the real number new_CSR format.

[0091] Figure 3 It is a schematic diagram for storing a pattern matrix in the complex CSR format.

[0092] Figure 4 It is a schematic diagram for splitting the Values array in the complex CSR format into the ValuesR array and the ValuesI array.

[0093] Figure 5 It is a schematic diagram for vectorizing the real number new_CSR format of 2 pattern matrices to obtain the fused new_CSR format.

[0094] Figure 6 It is an overall flowchart of an improved iterative solution algorithm. Detailed Embodiments

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

[0096] In the field of numerical analysis (such as the field of numerical analysis of electromagnetic fields, the field of numerical analysis of mechanics, etc.), there is often a need to solve matrices. Among these matrices to be solved (such as large dense complex matrices, large sparse complex matrices, etc.), there are many pattern matrices. Here, the pattern matrix refers to a matrix with the same non-zero element distribution pattern. For example, in the field of electromagnetic calculation, when the model and mesh discretization remain unchanged, the positions of non-zero elements in the equations obtained at different frequencies are the same. If these matrices are to be solved iteratively, each solution requires a call to the solution function, and it is impossible to perform parallel solution of two or more linear equations.

[0097] In this embodiment, for such large complex matrices with the same non-zero element distribution pattern (in fact, it can be applied not only to large complex matrices, but also to the iterative solution of all such pattern matrices. Generally, the direct solution of small matrices is more efficient), a data processing method for parallel iterative solution oriented to SIMD is designed. By using the SIMD extension component, simultaneous solution of multiple equations with the same non-zero element distribution in the horizontal direction is realized. For example, in the electromagnetic field, for a model that has been discretized into a grid, when solving the electromagnetic responses at two frequencies of 100 Hz and 10 Hz, a system of equations needs to be solved for each frequency, and the positions of non-zero elements in these two systems of equations are exactly the same, generating two pattern matrices. Then, it is suitable to use the method proposed by the present invention to vectorize the solution process in the horizontal direction to achieve parallel solution.

[0098] Please refer to Figure 1 , Figure 1 which is a flowchart of the data processing method for parallel iterative solution oriented to SIMD provided in this embodiment. In this embodiment, the data processing method for parallel iterative solution oriented to SIMD may include step S10, step S20, step S30, and step S40.

[0099] First, step S10 can be executed.

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

[0101] In this embodiment, these data to be processed, that is, the pattern matrices to be solved (usually including boundary conditions, written in the form of Ax = b), can be obtained. The number of pattern matrices processed at one time needs to satisfy 2 m , m ∈ Z +, each pattern matrix has the same non-zero element distribution pattern and is a complex matrix. Currently, the latest ARM9 platform supports the simultaneous execution of 4 double-precision instructions, and it may be expanded in the future. Although the technical difficulty is high and the number of double-precision instructions that can be processed synchronously is limited (currently at the level of simultaneously executing 2-4 double-precision instructions), there is still great application potential. Moreover, the idea of this solution can theoretically achieve a higher level of parallel processing to greatly improve the processing efficiency. For the convenience of description in this embodiment, two pattern matrices are used as examples for introduction, which should not be regarded as a limitation to this application.

[0102] After obtaining the data to be processed, step S20 can be executed.

[0103] Step S20: Store each pattern matrix in the real new_CSR format. Among them, the real 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 in each row of 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 the real new_CSR format, and the real new_CSR format includes the ValuesR array, the ValuesI array, the rowPtr array, and the colIndices array. As Figure 2 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 in each row of 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 there is already a CSR storage format in the prior art (for complex matrices, it is stored in the complex CSR format), the complex CSR format can be further processed and converted into the real new_CSR format.

[0106] As Figure 3As shown, exemplarily, each pattern matrix can be stored in a complex CSR format (which can be achieved by existing technologies). 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 in each row of the pattern matrix in the values array. The colIndices array is used to store the column indices of each non-zero element in the pattern matrix.

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

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

[0109] Get 2 m After obtaining the real new_CSR format corresponding to the two pattern matrices, step S30 can be executed.

[0110] Step S30: Vectorize the real new_CSR format of the two m pattern matrices to obtain a vectorized fused new_CSR format. Among them, the fused new_CSR format includes a new_ValuesR array, a new_ValuesI array, a rowPtr array, and a colIndices array. The new_ValuesR array includes several real part arrays, and each real part array is used to store the real parts of the non-zero elements at the same position of the two m pattern matrices. The new_ValuesI array includes the corresponding number of imaginary part arrays, and each imaginary part array is used to store the imaginary parts of the non-zero elements at the same position of the two m pattern matrices.

[0111] In this embodiment, the real new_CSR format of the two m pattern matrices can be vectorized. Specifically, the real parts of the non-zero elements at the same position in the ValuesR array of the real new_CSR format of the two m pattern matrices can be combined in the order of the pattern matrices to form a combined array containing 2 mThe real part arrays of the non-zero elements' real parts of the 2m pattern matrices are combined to form a new_ValuesR array containing several real part arrays; the imaginary parts of the non-zero elements at the same positions in the ValuesI array in the real new_CSR format of the 2m pattern matrices are combined in the order of the pattern matrices to form an imaginary part array containing 2 m non-zero element imaginary parts, obtaining a new_ValuesR array containing several imaginary part arrays; the rowPtr array remains unchanged; the colIndices array remains unchanged.

[0112] As Figure 5 shown, taking the vectorization of the real new_CSR format corresponding to 2 pattern matrices as an example, after converting the 2 pattern matrices into the corresponding real new_CSR format, vectorization is performed to obtain the corresponding fused new_CSR format.

[0113] After obtaining the fused new_CSR format corresponding to 2 m pattern matrices, 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] To implement the parallel iterative solution for 2 m pattern matrices, this embodiment designs an improved iterative solution algorithm, and the overall process is as Figure 6 shown.

[0116] First, perform initialization:

[0117] Use the SIMD extension component to read the quadruple in the fused new_CSR format of 2 m pattern matrices, the right values Breal and Bimag in the fused new_CSR format, the dimension n of the pattern matrices, the maximum number of iterations IMAX, and the convergence tolerance TOL. Among them, the quadruple in the fused new_CSR format is obtained by vectorizing the matrix A in the 2 m pattern matrices Ax = b stored in the real new_CSR format, and it contains the new_ValuesR array, the new_ValuesI array, the rowPtr array, and the colIndices array. The right value in the fused new_CSR format is obtained by vectorizing the right-hand vector b in the 2 m pattern matrices Ax = b stored in the real new_CSR format, and it contains 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, the search direction vectors Preal and Pimag, the residual vectors Rreal and Rimag, the solution vectors XR and XI, the residual norm VNRM, the errors RSS1, …, RSSi, …, RSS2 m , and assign values: Rreal = Breal, Rimag = Bima, Preal = Preal, Piamg = Rimage, 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 mode 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, 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 iteratively solved solution vectors XR and XI. Then decompose the solution vectors XR and XI (decompose them into solutions corresponding to 2 m mode matrices, and correspondingly merge the real and imaginary parts of the solutions of each mode matrix), and then the solution vectors of 2 m mode matrices can be obtained.

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

[0121] S1: Determine whether the current iteration termination condition is satisfied.

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

[0123] S2: If the iteration termination condition is satisfied, end the iteration and output the iteratively solved solution vectors XR and XI.

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

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

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

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

[0128] S34: Update the solution vectors XR and XI at the current iteration count based on the update step size α and the search direction vectors Preal and Pimag at the current iteration count.

[0129] S35: Update the residual vectors Rreal and Rimag based on the update step size α and the intermediate vectors AXR and AXI, to obtain the residual vectors Rreal′ and Rimag′ at the current iteration count.

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

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

[0132] S38: Increment the iteration count by 1, calculate the residual norm VNRM2 after iteration, and calculate the errors RSS1,…,RSSi…,RSS2 at the current iteration count based on the initialized residual norm VNRM and the residual norm VNRM2 after iteration. m , and then return to S1.

[0133] To introduce this solution in detail, the iteration process at this time (i.e., the iteration process at the current iteration count) is introduced in detail here.

[0134] First, in S31, it is necessary to call the improved cnorm function to record the inner product value RNORM of the residual vector at the current iteration count.

[0135] Initialize the vectors sum1 and sum2.

[0136] The SIMD extension component reads 2 m data from the residual vector Rreal, and reads 2 m data from the residual vector Rimag, where the step size for the SIMD extension component to read data is 2 m , Rreal is the real part of the residual vector R before update at the current iteration count, and Rimag is the imaginary part of the residual vector R before update at the current iteration count.

[0137] Then, calculate 2 according to the following formulam Inner product value RNORM of a residual vector:

[0138] RNORM = RTR, (1)

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

[0140] RNORMreal = RrealTRreal - RimagTRimag, (2)

[0141] RNORMimag = RrealTRimag + RimagTRreal, (3)

[0142] Where RNORMreal is the real part of the inner product value RNORM of the residual vector, RNORMimag is the imaginary part of the inner product value RNORM of the residual vector, RrealT is the transpose of Rreal, and RimagT is the transpose of Rimag.

[0143] After that, through CCR[0] = sum1[0], CCR[1] = sum1[1], …, CCR[k] = sum1[k], …, CCR[2 m -1] = sum1[2 m -1], the real part RNORMreal of the inner product value of 2 m residual vectors is saved. Through 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], the imaginary part RNORMimag of the inner product value of 2 m residual vectors is saved.

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

[0145] Taking the case of 2 mode matrices as an example, the improved cnorm function (i.e., the inner product function of the improved vector) is as follows:

[0146]

[0147] As shown above, the vdupq_n_f64, vld1q_f64, and vmlaq_f64 functions that use NEON registers are used for parallel computing to implement the dot product calculation of two sets of vectors at once. Thus, the calculation of the inner product value RNORM of the residual vector at the current iteration can be completed.

[0148] Secondly, in S32, the improved productAll function is called to calculate the intermediate vectors AXR and AXI at the current iteration.

[0149] The SIMD extension unit reads the search direction vectors Preal and Pimag, and the quadruples in the fused new_CSR format: the new_ValuesR array, the new_ValuesI array, the rowPtr array, and the colIndices array. Among them, the step size for the SIMD extension unit to read data is 2. m 。

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

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

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

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

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

[0155] 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, the rowPtr array, and the colIndices array through the SIMD extension unit. 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, the rowPtr array, and the colIndices array through the SIMD extension unit. Pimag is the imaginary part of the search direction vector P.

[0156] Taking the case of iterative solution of two pattern matrices by NEON as an example, the improved productAll function (i.e., the improved matrix-vector multiplication function) is as follows:

[0157]

[0158]

[0159] When implementing the steps of matrix-vector multiplication, only multiplying at the positions of non-zero values in the matrix-vector multiplication can save a lot of time. For example, the first for loop here is used to control the i-th row of the matrix, that is, select the i-th row to multiply the vector; the second for loop is used to control the multiplication calculation between 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. Thus, the calculation of the intermediate vectors AXR and AXI for the current iteration can be completed.

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

[0161] The SIMD extension unit reads the inner product value RNORMreal and RNORMimag of the residual vector, the search direction vectors Peak and Pimag, and the intermediate vectors AXR and AXI for the current iteration. Among them, the step size for the SIMD extension unit to read data is 2 m 。

[0162] Then, the update step size α for the current iteration can be calculated according to the following formula:

[0163]

[0164] Furthermore:

[0165]

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

[0167] Here, an improved complex multiplication function and an improved complex division function need to be used to implement.

[0168] Taking the case of iterative solution of 2 pattern matrices by NEON as an example, the improved complex multiplication function (ComMul_NEON function) is as follows:

[0169]

[0170] In the multiplication of complex numbers, NEON is used here to calculate the multiplication of two groups of complex numbers simultaneously to achieve parallel calculation. The input here is also virtual, because the complex multiplication will be used more than once in the whole iterative process. Each time the corresponding input is made, the parallel multiplication of complex numbers can be realized. The complex multiplication formula can be summarized as: for complex numbers Z1=(a + bi), Z2=(c + di), there is:

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

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

[0173]

[0174]

[0175] In the division of complex numbers, NEON is used here to calculate the division of two sets of complex numbers simultaneously to achieve parallel computing. The input here is also virtual because the complex division will be used in more than one place in the entire iterative process. Each time the corresponding input is given, the parallel division calculation of complex numbers can be realized. The complex division formula can be summarized as: for complex numbers Z1=(a + bi) and Z2=(c + di), there is:

[0176]

[0177] Thus, the calculation of the update step α at the current iteration can be completed. Note that the update step α here and the step length for the SIMD extension component to read data are 2 m are two different concepts and should not be confused.

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

[0179] The SIMD extension component reads the update step α at the current iteration and the search direction vectors Preal and Pimag. Among them, the step length for the SIMD extension component to read data is 2 m .

[0180] Then, calculate and update the solution vector at the current iteration according to the following formula:

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

[0182] According to the improved complex multiplication function, there is:

[0183] XR' = XR + αr * Preal - αi * Pimag, (13)

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

[0185] Where XR′ is the real part of the updated solution vector X at the current iteration, XR is the real part of the solution vector X before update at the current iteration, that is, the real part of the solution vector X updated in the previous 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 updated solution vector X at the current iteration, XI is the imaginary part of the solution vector X before update at the current iteration, that is, the imaginary part of the solution vector X updated in the previous iteration, αi is the imaginary part of the update step α, and Pimag is the imaginary part of the search direction vector P.

[0186] When updating the solution vector at the current iteration, the improved complex multiplication function (ComMul_NEON function) is used. For details, refer to the previous introduction and will not be elaborated here. Thus, the update of the solution vector at the current iteration can be realized.

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

[0188] Accordingly, the SIMD extension unit reads the update step α and the intermediate vectors AXR and AXI at the current iteration, where the step size for the SIMD extension unit to read data is 2 m .

[0189] Then, the residual vector at the current iteration can be updated according to the following formula:

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

[0191] According to the improved complex multiplication function, there is:

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

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

[0194] Among them, Rreal′ is the real part of the updated residual vector R at the current iteration, Rreal is the real part of the residual vector R before update at the current iteration, 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, Rimag′ is the imaginary part of the updated residual vector R at the current iteration, Rimag is the imaginary part of the residual vector R before update at the current iteration, that is, the imaginary part of the residual vector R updated in the previous iteration, and αi is the imaginary part of the update step α, and AXI is the imaginary part of the intermediate vector AX at the current iteration.

[0195] Similarly, when updating the residual vector at the current iteration, the improved complex multiplication function (ComMul_NEON function) is also used. For details, please refer to the previous introduction and will not be elaborated here. Thus, the update of the residual vector at the current iteration can be realized.

[0196] In S36, the conjugate parameter γ can be calculated based on the inner product value RNORM of the residual vector at the current iteration, and the real and imaginary parts Rreal′ and Rimag′ of the residual vector at the current iteration.

[0197] The SIMD extension component reads the inner product value RNORM of the residual vector at the current iteration, and the real and imaginary parts Rreal′ and Rimag′ of the residual vector at the current iteration. The step size for the SIMD extension component to read data is 2 m 。

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

[0199]

[0200] 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, there are:

[0201]

[0202] Among them, Rreal′ is the real part of the updated residual vector R′ at the current iteration, Rreal′T is the transpose of Rreal′, Rimag′ is the imaginary part of the updated residual vector R′ at the current iteration, Rimag′ T is the transpose of Rimag′, Rreal is the real part of the residual vector R before update at the current iteration, that is, the real part of the residual vector R updated in the previous iteration, Rreal T is the transpose of Rreal, Rimag is the imaginary part of the residual vector R before update at the current iteration, that is, the imaginary part of the residual vector R updated in the previous iteration, Rimag TIs the transpose of Rimag.

[0203] When calculating the conjugate parameter γ at the current iteration, the improved complex multiplication function (ComMul_NEON function) and the improved complex division function (ComDivi_NEON function) are also used. For details, please refer to the previous introduction and will not be elaborated here. Thus, the calculation of the conjugate parameter γ at the current iteration can be realized.

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

[0205] 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 iteration, and the residual vectors Rreal′ and Rimag′ at the current iteration (note that they are the updated residual vectors at the current iteration). Among them, the step size for the SIMD extension component to read data is 2. m .

[0206] Then, update the search direction vectors Preal and Pimag according to the following formula:

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

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

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

[0210] Pimag′ = Rimag′ + γr * Pimag + yi * Preak, (23)

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

[0212] When updating the search direction vector at the current iteration, the improved complex multiplication function (ComMul_NEON function) is also used. For details, please refer to the previous introduction and will not be elaborated here. Thus, the update of the search direction vector at the current iteration can be realized.

[0213] In S38, after the above iterative update, the iteration count is incremented by 1, i.e., iter′ = iter + 1. Then, the residual norm VNRM2 after iteration can be calculated.

[0214] The SIMD extension unit reads the residual vectors Rreal′ and Rimag′ at the current iteration count, where the step size for the SIMD extension unit to read data is 2 m 。

[0215] Calculate the residual norm VNRM2 of each pattern matrix after this iteration according to the following formula:

[0216]

[0217] where VNRM2 is the residual norm calculated for the pattern matrix at the current iteration count, R′ is the updated residual vector at the current iteration count, n is the dimension of the number of rows or columns in the pattern matrix, Rreal′ j is the j-th value of Rreal′, Rimag′ j is the j-th value of Rimag′.

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

[0219]

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

[0221] After that, the errors RSS1,…,RSSi…,RSS2 at the current iteration count can be further calculated m 。

[0222] The SIMD extension unit reads the initialized residual norm VNRM and the residual norm VNRM2 after iteration, where the step size for the SIMD extension unit to read data is 2 m 。

[0223] Then, calculate the error of each pattern matrix at the current iteration count according to the following formula:

[0224]

[0225] where RSSi is the error of the i-th pattern matrix at the current iteration number, and i ∈ [1, 2 m .

[0226] The above are all the processes of an iteration cycle. After completing the iterative update of this cycle, it is possible to return to S1 to determine whether the iteration termination condition is satisfied. The iteration termination condition is designed as: the iteration number iter reaches the maximum iteration number IMAX, or the error RSS of each pattern matrix at the current iteration number is not greater than the convergence tolerance TOL.

[0227] To verify the effect of this solution, a test case is provided here (in this embodiment, the NEON register is taken as an example, as the SIMD extension component, to run the iterative process in the data processing method for parallel iterative solution oriented to SIMD):

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

[0229] Unit / Method NEON Parallel Iterative Solution Iterative Solution NEON Efficiency Improvement Time (s) 2696.53 8087.68 66.66％

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

[0231] In summary, the embodiment of this application provides a data processing method for parallel iterative solution oriented to SIMD, which uses the SIMD extension component to design a parallel iterative solution scheme for 2 m pattern matrices to achieve accelerated solution. Specifically, by obtaining the data to be processed (including 2 m pattern matrices, m ∈ Z + , each pattern matrix has the same non-zero element distribution pattern and is a complex matrix); storing each pattern matrix in the real new_CSR format, where the real 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 in each row of 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 mVectorize the real number new_CSR format of a pattern matrix to obtain a vectorized fused new_CSR format. The fused new_CSR format includes a new_ValuesR array, a new_ValuesI array, a rowPtr array, and a colIndices array. The new_ValuesR array contains several real part arrays, and each real part array is used to store the real parts of the non-zero elements at the same position of 2 m pattern matrices. The new_ValuesI array contains a corresponding number of imaginary part arrays, and each imaginary part array is used to store the imaginary parts of the non-zero elements at the same position of 2 m pattern matrices; Use the SIMD extension component to read the data in the vectorized fused new_CSR format and perform parallel iterative solution, and output the iterative solution result. This solution improves the storage format of CSR. By vectorizing and storing the corresponding data of multiple pattern matrices, the SIMD extension component (such as the NEON register) can read data in parallel and execute SIMD instructions (such as NEON instructions), and improve the relevant iterative solution algorithm, so as to achieve 2 m parallel iterative solution of pattern matrices, greatly improving the iterative solution efficiency.

[0232] By storing each pattern matrix in the complex CSR format, 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 in each row of 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 in the Values array of the complex CSR format into the real part of the non-zero element and the imaginary part of the non-zero element, store the real parts of the non-zero elements in order to obtain the new_ValuesR array, store the imaginary parts of the non-zero elements in order to obtain the new_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 processing to obtain a new real number new_CSR format, which can make full use of the existing technology and be applied to a wide range of application scenarios. And the 2 m real parts of the non-zero elements at the same position in the new_ValuesR array of the real number new_CSR format of pattern matrices are combined into a real part array containing 2 m real parts of non-zero elements to obtain a new_ValuesR array containing several real part arrays; The 2 mIn the real number new_CSR format of a pattern matrix, the imaginary parts of the non-zero elements at the same position in the ValuesI array are combined in the order of the pattern matrix to form an imaginary part array containing 2 m imaginary part arrays of non-zero element imaginary parts, resulting in a new_ValuesR array containing several imaginary part arrays; the rowPtr array remains unchanged; the colIndices array remains unchanged. Since accelerating the program using vector registers requires reading the data in memory into them in advance, vectorization is performed based on the real number new_CSR format to obtain a fused new_CSR format, which can achieve vectorization and make the fused new_CSR format adapt to the SIMD extension components. The alignment of data can significantly improve the performance of the SIMD extension components, realize parallel iterative solution calculations, effectively improve the calculation speed of parallel processing, greatly improve the calculation efficiency, and shorten the iterative solution time.

[0233] The parallel iterative solution scheme provided in this embodiment improves the iterative solution algorithm, enabling the entire iterative solution algorithm to be applicable to the calculation of parallel iterative solutions for 2 m pattern matrices. Corresponding improved formulas are designed to achieve parallel iterative updates. At the same time, the operation step (such as reading data) of the SIMD extension components is adjusted to ensure the stable operation of the parallel iterative solution. The test case uses this scheme (taking parallel acceleration iterative solution with NEON registers as an example) to perform iterative solutions on two complex sparse matrices of size 122036 * 122036. The non-zero value of each matrix is 8168600. The comparison method is iterative solution without using SIMD extension components for acceleration. CG iterative solutions are performed on the two matrices respectively, and the obtained solution times are 2696.53 seconds and 8087.68 seconds respectively. The iterative solution efficiency is improved by 66.66%. If more matrices are solved, the acceleration performance will be more obvious.

[0234] The above are only the embodiments of the present application and are not used to limit the protection scope of the present application. For those skilled in the art, various changes and modifications can be made to the present application. 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 parallel iterative solution oriented to SIMD, characterized in that Including: Obtain the data to be processed, where the data to be processed contains 2 m pattern matrices, m ∈ Z + , and each pattern matrix has the same non-zero element distribution pattern and is a complex matrix; Storing each pattern matrix in the real new_CSR format, where the real new_CSR format includes a ValuesR array, a ValuesI array, a rowPtr array, and a colIndices array. The ValuesR array is used to store the real parts of the non-zero elements of the pattern matrix, the ValuesI array is used to store the imaginary parts of the non-zero elements of the pattern matrix, the rowPtr array is used to store the positions of the first non-zero element in each row of the pattern matrix in the values array, and the colIndices array is used to store the column indices of each non-zero element in the pattern matrix; For 2 m Vectorize the real new_CSR format of two pattern matrices to obtain a vectorized fused new_CSR format. The fused new_CSR format includes a new_ValuesR array, a new_ValuesI array, a rowPtr array, and a colIndices array. The new_ValuesR array contains a number of real part arrays, and each real part array is used to store the real parts of the non-zero elements at the same position of two m pattern matrices. The new_ValuesI array contains a corresponding number of imaginary part arrays, and each imaginary part array is used to store the imaginary parts of the non-zero elements at the same position of two m pattern matrices; Using the SIMD extension component to read the vectorized fused new_CSR format data, performing parallel iterative solution, and outputting the iterative solution result; Among them, using the SIMD extension component to read the vectorized fused new_CSR format data and performing parallel iterative solution includes: Read using SIMD extension components 2 m Fused quadruples in new_CSR format of 2 pattern matrices, right values Breal and Bimag in fused new_CSR format, dimension n of the pattern matrix, maximum number of iterations IMAX, convergence tolerance TOL, where the fused quadruples in new_CSR format are 2 m Obtained by vectorizing the matrix A in the 2 pattern matrices Ax = b after being stored in real new_CSR format, including new_ValuesR array, new_ValuesI array, rowPtr array, and colIndices array, and the right value in fused new_CSR format is 2 m Obtained by vectorizing the right-hand vector b in the 2 pattern matrices Ax = b after being stored in real new_CSR format, including Breal and Bimag, where Breal is the real part of the right value and Bimag is the imaginary part of the right value; Initialize the intermediate vectors AXR and AXI, the search direction vectors Preal and Pimag, the residual vectors Rreal and Rimag, the solution vectors XR and XI, the residual norm VNRM, and the errors RSS1, …, RSSi, …, RSS2 m , and assign values: Rreal = Breal, 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 respectively, Preal and Pimag are the real and imaginary parts of the search direction vector respectively, Rreal and Rimag are the real and imaginary parts of the residual vector respectively, and RSSi is the error corresponding to the i-th mode 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, and the errors RSS1, …, RSSi, …, RSS2 m , the solution vectors XR and XI are iteratively updated using the SIMD extension component, and finally the iteratively solved solution vectors XR and XI are output; Decompose the solution vectors XR and XI to obtain the solution vectors of m two pattern matrices.

2. The data processing method for SIMD-oriented parallel iterative solution according to claim 1, characterized in that Storing each pattern matrix in the real new_CSR format, including: Storing each pattern matrix in the complex CSR format, where 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 positions of the first non-zero element in each row of the pattern matrix in the values array, and the colIndices array is used to store the column indices of each non-zero element in the pattern matrix; For the complex CSR format of each pattern matrix: splitting the complex non-zero elements in the Values array of the complex CSR format into the real part of the non-zero element and the imaginary part of the non-zero element, storing the real part of the non-zero element in order to obtain the ValuesR array, storing the imaginary part of the non-zero element in order to obtain the ValuesI array, and keeping the rowPtr array and the colIndices array unchanged.

3. The data processing method for SIMD-oriented parallel iterative solution according to claim 1, wherein Pair 2 m Vectorize the real number new_CSR format of the 2 pattern matrices to obtain the vectorized fused new_CSR format, including: Take 2 m The real parts of the non-zero elements at the same position in the ValuesR array in the real new_CSR format of 2 pattern matrices are combined in the order of the pattern matrices into a real part array containing 2 m non-zero element real parts, obtaining a new_ValuesR array containing several real part arrays; Take 2 m For the non-zero element imaginary parts at the same position in the ValuesI array of the real new_CSR format of 2 pattern matrices, they are combined in the order of the pattern matrices into an imaginary part array containing 2 m non-zero element imaginary parts, resulting in a new_ValuesR array containing several imaginary part arrays; The rowPtr array remains unchanged; The colIndices array remains unchanged.

4. The data processing method for SIMD-oriented parallel iterative solution according to claim 1, wherein Based on the intermediate vectors AXR and AXI, search direction vectors Preal and Pimag, residual vectors Rreal and Rimag, residual norm VNRM, errors RSS1, …, RSSi, …, RSS2 m , the solution vectors XR and XI are iteratively updated using the SIMD extension component, and finally the iteratively solved solution vectors XR and XI are output, including: S1: Judging whether the current iteration termination condition is satisfied; S2: If the iteration termination condition is satisfied, end the iteration and output the solution vectors XR and XI of the iterative solution; S3: If the iteration termination condition is not satisfied, start the current iteration: S31: Calling the improved cnorm function to record the inner product value RNORM of the residual vector at the current iteration; S32: Calling the improved productAll function to calculate the intermediate vectors AXR and AXI at the current iteration; S33: Calculating the update step size α at the current iteration based on the inner product value RNORM of the residual vector and the intermediate vectors AXR and AXI at the current iteration; S34: Updating the solution vectors XR and XI at the current iteration based on the update step size α at the current iteration and the search direction vectors Preal and Pimag; S35: Update the residual vectors Rreal and Rimag based on the update step α, the intermediate vectors AXR and AXI at the current iteration number to obtain the residual vectors Rreal' and Rimag' at the current iteration number; S36: Calculate the conjugate parameter γ based on the inner product value RNORM of the residual vectors 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: Increment the iteration count by 1, calculate the residual norm VNRM2 after iteration, and calculate the errors RSS1, …, RSSi, …, RSS2 at the current iteration count based on the initialized residual norm VNRM and the residual norm VNRM2 after iteration, m and then return to S1.

5. The data processing method for SIMD-oriented parallel iterative solution according to claim 4, characterized in that Call the improved cnorm function and record the inner product value RNORM of the residual vectors at the current iteration number, including: Initialize the vectors sum1 and sum2; The SIMD extension unit reads 2 m data from the residual vector Rreal and reads 2 m data from the residual vector Rimag. Here, the step size for the SIMD extension unit to read data is 2 m , where Rreal is the real part of the residual vector R before update at the current iteration number, and Rimag is the imaginary part of the residual vector R before update at the current iteration number; Calculate according to the following formula 2 m Inner product value RNORM of residual vectors: RNORM = R T R, RNORMreal = Rreal T Rreal - Rimag T Rimag, RNORMimag = Rreal T Rimag + Rimag T Rreal, Where RNORMreal is the real part of the inner product value RNORM of the residual vectors, RNORMimag is the imaginary part of the inner product value RNORM of the residual vectors, 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] to save the real part RNORMreal of the inner product values of 2 m residual vectors. 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] to save the imaginary part RNORMimag of the inner product values of 2 m residual vectors.

6. The data processing method for SIMD-oriented parallel iterative solution according to claim 4, characterized in that, Call the improved productAll function to calculate the intermediate vectors AXR and AXI at the current iteration number, including: The SIMD extension unit reads the search direction vectors Preal and Pimag, and the quadruples in the fused new_CSR format: the new_ValuesR array, the new_ValuesI array, the rowPtr array, and the colIndices array. Here, the step size for the SIMD extension unit to read data is 2 m ; Calculate the intermediate vectors AXR and AXI at the current iteration number according to the following formula: AXR = AR * Preal - AI * Pimag, AXI = AR * Pimag + AI * Preal, Where AXR is the real part of the intermediate vector AX, AR is the real part of the pattern matrix, 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, 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.

7. The data processing method for SIMD-oriented parallel iterative solution according to claim 5, characterized in that Calculate the update step α at the current iteration number based on the inner product value RNORM of the residual vectors at the current iteration number and the intermediate vectors AXR and AXI, including: The SIMD extension unit reads the inner product values RNORMreal and RNORMimag of the residual vector, the search direction vectors Preal and Pimag, and the intermediate vectors AXR and AXI at the current iteration number. Herein, the step size for the SIMD extension unit to read data is 2 m ; Calculate the update step α at the current iteration number according to the following formula: Furthermore: Where αr is the real part of the update step α, and αi is the imaginary part of the update step α.

8. The data processing method for SIMD-oriented parallel iterative solution according to claim 7, wherein Update the solution vectors XR and XI at the current iteration number based on the update step α at the current iteration number and the search direction vectors Preal and Pimag, including: The SIMD extension unit reads the update step α, the search direction vectors Preal and Pimag at the current iteration count. Herein, the step size for the SIMD extension unit to read data is 2 m ; Calculate the updated solution vectors at the current iteration number 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 updated solution vector X at the current iteration, XR is the real part of the solution vector X before update at the current iteration, that is, the real part of the solution vector X updated in the previous iteration, αr is the real part of the update step size α, Preal is the real part of the search direction vector P, XI′ is the imaginary part of the updated solution vector X at the current iteration, XI is the imaginary part of the solution vector X before update at the current iteration, that is, the imaginary part of the solution vector X updated in the previous iteration, and αi is the imaginary part of the update step size α, Pimag is the imaginary part of the search direction vector P.

9. The data processing method for SIMD-oriented parallel iterative solution according to claim 7, wherein Based on the update step size α, and the intermediate vectors AXR and AXI at the current iteration, update the residual vectors Rreal and Rimag to obtain the residual vectors Rreal′ and Rimag′ at the current iteration, including: The SIMD extension component reads the update step α, the intermediate vectors AXR and AXI at the current iteration count. Here, the step size for the SIMD extension component to read data is 2 m ; Calculate and update the residual vectors at the current iteration according to the following formula: Rreal′ = Rreal - αr * AXR + αi * AXI, Rimag′ = Rimag - αr * AXI - αi * AXR, where, Rreal′ is the real part of the residual vector R updated at the current iteration, Rreal is the real part of the residual vector R before update at the current iteration, that is, the real part of the residual vector R updated in the previous iteration, αr is the real part of the update step size α, AXR is the real part of the intermediate vector AX at the current iteration, Rimag′ is the imaginary part of the residual vector R updated at the current iteration, Rimag is the imaginary part of the residual vector R before update at the current iteration, that is, the imaginary part of the residual vector R updated in the previous iteration, αi is the imaginary part of the update step size α, and AXI is the imaginary part of the intermediate vector AX at the current iteration.

Images