Vectorization parallel solving method and device for frequency domain electromagnetic forward modeling calculation

Abstract

The invention provides a vectorization parallel solving method and device for frequency domain electromagnetic forward modeling calculation, and relates to the technical field of algorithm solving. A large-scale sparse linear system obtained after discretization in the frequency domain electromagnetic forward modeling calculation process is converted into a vectorization storage format; sparse linear systems of multiple frequencies are solved at the same time through a vector instruction of a modern CPU by utilizing a batch processing thought. When the sparse linear system is solved by using the generalized minimum residual error method, expanding each step of core operator into frequency-dependent vector calculation so as to achieve the target of solving the sparse linear system with multiple frequencies at one time; according to the method, the calculation characteristics and the memory access characteristics of the multi-frequency sparse linear system are combined, the long vector register of the CPU is fully utilized, the memory access frequency of the matrix is reduced, the problems of low calculation efficiency and low calculation resource utilization rate in the prior art are solved, quick and efficient multi-frequency frequency domain electromagnetic forward modeling calculation is realized, and the calculation efficiency is greatly improved.

Description

Technical Field

[0001] The present invention relates to the technical field of algorithm solving, and particularly to a vectorized parallel solving method and device for frequency-domain electromagnetic forward calculation. Background Art

[0002] Frequency-domain electromagnetic forward calculation is one of the core technologies in geophysical exploration (especially electromagnetic methods); frequency-domain electromagnetic forward calculation can simulate the propagation behavior of electromagnetic waves in underground media at different frequencies for exploration objects, so as to obtain key information such as resistivity distribution, electrical layer thickness, geometric shape, electrical interface position, phase information, and the size and position of abnormal bodies of the underground structure, thereby obtaining the electromagnetic response characteristics of the exploration object. These electromagnetic response characteristics are of great significance for fields such as geological exploration, resource investigation, and environmental monitoring.

[0003] Generally, after discretization by the finite element method (FEM) or the finite difference time domain method (FDTD) in forward calculation, a large-scale sparse linear system is usually generated. The large-scale sparse linear system usually occupies the main calculation time of the forward calculation, usually more than 80%. The generalized minimum residual method (GMRES) is an iterative algorithm for efficiently solving large-scale sparse non-symmetric linear systems and is widely used in electromagnetic field forward calculation. Computational performance and time consumption are crucial for actual production efficiency. Therefore, it is necessary to improve the forward calculation speed by combining the characteristics of the modern general-purpose processor CPU architecture.

[0004] Currently, the commonly used calculation methods for frequency-domain electromagnetic forward calculation are one is to solve the linear systems generated by frequencies sequentially, and the other is to solve them in parallel by allocating different frequencies to multiple CPU cores.

[0005] Currently, CPU processors usually provide vector registers of 256 - 512 bits, which can calculate 4 - 8 64-bit floating-point calculations or 8 - 16 32-bit floating-point calculations simultaneously. However, the currently commonly used calculation methods do not combine the calculation characteristics and memory access characteristics of multi-frequency sparse linear systems, not only will repeatedly read the sparse linear matrix, but also cannot give full play to the characteristics of the long vectorization of the current CPU processors, resulting in the problems of reduced calculation efficiency and low utilization rate of calculation resources.

[0006] Therefore, how to efficiently utilize the long vector registers of the CPU and how to improve the memory access efficiency of the matrix to solve multi-frequency frequency-domain electromagnetic forward calculation is a difficult problem to be solved urgently. Summary of the Invention

[0007] In order to overcome the deficiencies in the background art, the present invention provides a vectorized parallel solving method and device for frequency-domain electromagnetic forward calculation.

[0008] In order to achieve the above-mentioned invention object, the present invention adopts the following technical scheme: In a first aspect, the present invention provides a vectorized parallel solution method for frequency domain electromagnetic forward modeling, comprising the following steps: S1. According to the vector register length and calculation accuracy of the CPU processor, the sparse linear system A with different frequencies is i x i = b i Convert to vectorized computing format, including sparse matrix A i The real matrix Re and the sparse matrix A i The imaginary matrix Im i , reference complex matrix Im, frequency value f i , frequency vector f vec and right end b i_vec Among them, A i is a sparse matrix of different frequencies, A i for nrow × nrow The complex matrix of ; x i is the vector to be solved at different frequencies, with dimensions of nrow * 1 vector of; b i is the right-hand term vector of different frequencies, with dimension nrow * 1 Vector ；i are the serial numbers of different frequencies; i =1,2,……n, n is the number of frequencies, n is a positive integer greater than or equal to 2; nrow is the number of degrees of freedom that need to be solved; S2. Initialization x i The approximate solution of x i_vec , and use vectorized calculation method to calculate the residual r i_vec and its bi-norm β vec Post-convergence update of approximate solution x i_vec ; S3. Initialize the orthogonal basis matrix V vec and the Hessenberg matrix H vec , and calculate V vec [:, 0], V vec [:,0] = ri_vec / β vec ; S4. Update the orthogonal basis matrix V vec and the Hessenberg matrix H vec by iterative orthogonalization, Hessenberg matrix update, and least squares solution, and calculate the current residual Res vec . When the current residual Res vec converges or the maximum number of iterations is reached, calculate the updated approximate solution V vec according to the orthogonal basis matrix x i_vec .

[0009] Specifically, step S1 specifically includes the following steps: S11. Obtain the vector register length of the CPU processor; S12. Calculate the number of vectors that can be calculated simultaneously by one machine instruction of the CPU according to the vector register length and the calculation accuracy; S13. Convert the sparse matrix A i with different frequencies into A i =(Re + Im i ), and then convert Im i into Im i = Im· f i ; where Re is the real matrix of the sparse matrix A i , and Im i is the imaginary matrix of the sparse matrix A i ; Im represents the reference complex matrix, f i is the frequency value; S14. Construct the frequency vector f i through the frequency value f vec , and pad f vec to a multiple of the number of vectors. After padding, the vector length of the frequency vector f vec is vec_num; S15. Construct the right-hand side b i through the right-hand side vector b i_vec , b i_vec is a vector of nrow rows, and pad b i_vecEach row is padded to a multiple of the number of vectors, and the padded value is the right-hand side vector b i Any element in

[0010] Specifically, step S2 includes the following steps: S21. Set the maximum number of iterations max_iter and the convergence tolerance tolerance ; S22. Initialize x i approximate solution of x i_vec , x i_vec which is a nrow vector with vec_num rows, and each element in it is a vector with r i_vec length; r i_vec = b i_vec -(Re + Im· f i_vec ) x i_vec ; r i_vec is a nrow vector with vec_num rows, and each element in it is a S24. Calculate the initial residual r i_vec two-norm of β vec : β vec = sqrt ( r 0_vec 2 + r 1_vec 2 + r 2_vec 2 +…… r n_vec 2 ) where r n_vec represents a set of vectors, β vec is a vector with dimension vec_num , and n = nrow - 1

[0011] Specifically, step S3 specifically includes the following steps: S31. Initialize the orthogonal basis matrix V vec , where the orthogonal basis matrix V vec is a matrix of dimension nrow ×( max_iter +1), and each matrix element is a vector of length vec_num ; S32. Calculate V vec [:, 0], V vec [:, 0] = r i_vec / β vec ; S33. Initialize the Hessenberg matrix H vec , where the matrix H vec is constructed as a matrix of dimension max_iter ×( max_iter +1), and each matrix element is a vector of length vec_num ;

[0012] Specifically, step S4 specifically includes the following steps: S41. Set the current iteration loop count iter = 0; S42. Calculate the orthogonal vector w i_vec = (Re + Im· f i_vec ) * V vec [:, iter , w i_vec is a vector of nrow rows, and each element therein is a vector of length vec_num ; S43. Perform modified orthogonalization; S44. Calculate the two-norm of the orthogonal vector w i_vec and assign it to the Hessenberg matrix H vec : H vec iter +1] iter = sqrt ( w 0_vec 2 + w 1_vec ​​2 + w 2_vec 2 +…… w n_vec 2 )), where n = nrow -1; S45. Standardize the Hessenberg matrix H vec as the orthogonal basis matrix V vec : V vec [:, iter + 1] = w i_vec / H vec iter +1] iter ; S46. Solve the Hessenberg matrix H vec by the least squares method to obtain the first vector array y i_vec and the second vector array e i_vec ; S47. Calculate the current residual Res vec based on the Hessenberg matrix H y i_vec the first vector array e i_vec and the second vector array vec and check for convergence; S48. If the current residual Res vec has not converged, then determine whether iter is less than the maximum number of iterations max_iter. If so, jump to step S42; otherwise, jump to step S49; S49. Update the approximate solution V vec based on the orthogonal basis matrix y i_vec and the first vector array x i_vec to obtain the solution of the frequency-domain electromagnetic forward modeling: x i_vec = x i_vec + V vec [:, iter + 1] y i_vec .

[0013] Specifically, step S43 specifically includes the following steps: S431. Set​k = 0; Enter the loop; S432. Projection coefficient H vec [k] iter Calculate: H vec [k] iter = dot (( V vec [:, k ) T , w i_vec ), dot indicating the sum of the dot products of two vectors; S433. Through w i_vec = w i_vec - H vec k iter V vec [:, k for orthogonalization correction; S434. If k is less than iter + 1 ,k = k + 1, jump to S432, otherwise jump to S435; S435. The orthogonalization correction ends.

[0014] Specifically, step S46 specifically includes the following steps: S461. Set k = 0; Enter the loop; S462. Construct H k as a scalar matrix of ( iter + 2) × ( iter + 1), and store the k H vec [: iter + 2, : iter + 1] in the k +1-th element of each vector element; S463. Construct e k as a scalar array of ( iter + 2), initialize e k to 0; e k [0] = β k ; S464. Solve for H by the least squares method k y k = e​​​​​k , obtain y k which is a scalar array of ( iter +1); S465. If k is less than vec_num,k = k +1 ， jump to S462; otherwise, jump to S466; S466. Assemble y k (k = 0, 1, …… vec_num -1) into a new first vector array y i_vec , y i_vec which is a iter +1-row vector, and each element of which is a vec_num -long vector; S467. Assemble e k (k = 0, 1, …… vec_num -1) into a new second vector array e i_vec , e i_vec which is a iter +2-row vector, and each element of which is a vec_num -long vector.

[0015] Specifically, step S47 specifically includes the following steps: S471. Calculate tmp vec = H vec [: iter +2, : iter +1] y i_vec -e i_vec , tmp vec which is a iter +2-row vector, and each element of which is a vec_num -long vector; S472. Calculate Res vec , and the residual value of the current iteration is as follows. Res vec is a vector of vec_num length: Res vec = sqrt ( tmp 0_vec 2 + tmp 1_vec2 + tmp 2_vec 2 +…… tmp n_vec 2 ); S473. Determine whether all values in the current residual Res vec are less than the convergence tolerance tolerance . If so, jump to step S49; otherwise, jump to step S48.

[0016] Specifically, the sparse linear system A x = b is obtained by discretizing the frequency-domain electromagnetic forward calculation through the finite element method or the finite difference method.

[0017] Specifically, in step S4, it further includes: obtaining the electromagnetic response characteristics of the exploration object according to the solution of the frequency-domain electromagnetic forward calculation.

[0018] In a second aspect, the present invention provides a vectorized parallel solving device for frequency-domain electromagnetic forward calculation, including the following units: A vectorization unit, configured to convert sparse linear systems A at different frequencies i x i = b i into a vectorized calculation format, including the real matrix Re of the sparse matrix A i , the imaginary matrix Im of the sparse matrix A i , the reference complex matrix Im, the frequency value i , the frequency vector f i , and the right-hand side term f vec ; where A b i_vec is a sparse matrix at different frequencies, and A i is a complex square matrix of i × nrow × nrow ; x i is the vector to be solved at different frequencies, with a dimension of nrow * 1 vector; b i is the right-hand side term vector at different frequencies, with a dimension of nrow * 1 vector ；i is the serial number at different frequencies; i = 1, 2, …… n, where n is the number of frequencies, and n is a positive integer greater than or equal to 2; nrow is the number of degrees of freedom to be solved; The first initialization unit is used for initialization x i of the approximate solution x i_vec , and calculates the residual by using the vectorized calculation method r i_vec and its two-norm β vec and then converges to update the approximate solution x i_vec ; The second initialization unit is used for initializing the orthogonal basis matrix V vec and the Hessenberg matrix H vec , and calculates V vec [:, 0], V vec [:, 0] = r i_vec / β vec ; The iterative calculation unit is used for updating the orthogonal basis matrix V vec and the Hessenberg matrix H vec through iterative orthogonalization, Hessenberg matrix update and least squares solution, and calculating the current residual Res vec , and when the current residual Res vec converges or the number of iterations reaches the maximum, calculates and updates the approximate solution V vec according to the orthogonal basis matrix x i_vec , so as to obtain the solution of the frequency-domain electromagnetic forward calculation.

[0019] In a third aspect, the present invention also provides an electronic device, including a processor, a memory, a communication interface, and one or more programs, the one or more programs are stored in the memory and are configured to be executed by the processor, and the programs include instructions for executing the steps in any one of the methods in the first aspect.

[0020] The present invention proposes a vectorized parallel solution method for frequency-domain electromagnetic forward calculation, including: S1. Converting sparse linear systems A at different frequencies i x i = b i into a vectorized calculation format according to the vector register length and calculation accuracy of the CPU processor; S2. Initializing x i of the approximate solutionx i_vec , and calculate the residual using the vectorized calculation method r i_vec and its two-norm β vec and then converge to update the approximate solution x i_vec ; S3. Initialize the orthogonal basis matrix V vec and the Hessenberg matrix H vec , and calculate V vec [:, 0], V vec [:, 0] = r i_vec / β vec ; S4. Update the orthogonal basis matrix V vec and the Hessenberg matrix H vec by iterative orthogonalization, Hessenberg matrix update and least squares solution, and calculate the current residual Res vec , and when the current residual Res vec converges or the maximum number of iterations is reached, calculate and update the approximate solution V vec according to the orthogonal basis matrix x i_vec , so as to obtain the solution of the frequency-domain electromagnetic forward calculation, and obtain the electromagnetic response characteristics of the exploration object according to the solution of the frequency-domain electromagnetic forward calculation. The present invention converts the left-end matrix and the right-end term of the large-scale sparse linear system obtained after discretization in the frequency-domain electromagnetic forward calculation process into a vectorized storage format. Through the idea of batch processing, solve the sparse linear systems of multiple frequencies simultaneously, and calculate through the vector instructions of modern CPUs. When using the generalized minimum residual method to solve the sparse linear system, expand each step of the core operator therein (mainly including calculations such as sparse matrix-vector multiplication, vector dot multiplication, and vector scalar multiplication) into frequency-related vector calculations, and set the convergence criterion based on vectorized calculation to achieve the goal of solving the sparse linear systems of multiple frequencies at one time; the present invention combines the calculation characteristics and memory access characteristics of the multi-frequency sparse linear system, makes full use of the long vector registers of the CPU, reduces the memory access times of the matrix, solves the problems of reduced calculation efficiency and low utilization rate of calculation resources in the prior art, realizes fast and efficient multi-frequency frequency-domain electromagnetic forward calculation, and greatly improves the calculation efficiency of the frequency-domain electromagnetic forward calculation.

[0021] In addition, compared with traditional frequency parallel computing, this embodiment does not require any additional computing resources. The saved CPU cores can be used for other computations, such as parallel computing of matrix rows or computations of other programs, which can save more computing resources.

[0022] During metal ore exploration, metal ore bodies usually have significant electrical anomalies (such as low resistivity or high resistivity), but their burial depths, shapes, and scales are complex and variable. Traditional forward modeling methods are time-consuming in calculation and difficult to meet the needs of large-scale exploration. The vectorized parallel solution method for frequency-domain electromagnetic forward modeling of the present invention can greatly improve the data processing speed, quickly identify electromagnetic response characteristics such as the position of the ore body and electrical parameters, shorten the data processing cycle by more than 50%, significantly reduce the hardware cost, improve the accuracy of ore body boundary identification at the same time, reduce the blindness of drilling, and reduce the exploration cost.

[0023] During oil and gas exploration, oil and gas reservoirs usually show high-resistivity anomalies, but their electromagnetic responses are affected by complex formation structures (such as salt domes, faults). Traditional forward modeling methods are difficult to quickly simulate electromagnetic response characteristics such as the electromagnetic field distribution at multiple frequency points and multiple scales. The vectorized parallel solution method for frequency-domain electromagnetic forward modeling of the present invention can reduce the dependence on high-performance computing clusters, efficiently simulate electromagnetic responses under complex geological conditions, shorten the data processing cycle by more than 50%, significantly reduce the hardware cost, improve the accuracy of oil and gas reservoir boundary identification at the same time, and reduce the exploration risk. Brief Description of the Drawings

[0024] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required in the embodiments. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0025] Figure 1 is a schematic diagram of a vectorized parallel solution method for frequency-domain electromagnetic forward modeling provided according to an embodiment of the present invention; Figure 2 is a schematic diagram of a vectorized parallel solution device for frequency-domain electromagnetic forward modeling provided according to an embodiment of the present invention; Figure 3 is a schematic diagram of a vectorized parallel solution device for frequency-domain electromagnetic forward modeling provided according to an embodiment of the present invention. Detailed Embodiments

[0026] The present invention can be explained in detail by the following embodiments. The purpose of providing the present invention is to protect all technical improvements within the scope of the present invention. In the description of the present invention, it should be understood that if there are terms such as "upper", "lower", "front", "rear", "left", "right", etc. indicating the orientation or positional relationship, it is only corresponding to the drawings of the present application. For the convenience of describing the present invention, it does not indicate or imply that the device or element referred to must have a specific orientation.

[0027] Embodiment 1

[0028] Reference Figure 1 , this embodiment provides a vectorized parallel solution method for frequency-domain electromagnetic forward calculation, including the following steps: S1. According to the vector register length and calculation accuracy of the CPU processor, convert the sparse linear system A at different frequencies i x i = b i into a vectorized calculation format, including the real matrix Re of the sparse matrix A i , the imaginary matrix Im of the sparse matrix A i , the reference complex matrix Im, the frequency value i , the frequency vector f i , and the right-hand side term f vec ; where A b i_vec is a sparse matrix at different frequencies, i is the vector to be solved at different frequencies, x i is the right-hand side vector at different frequencies, b i is the serial number at different frequencies i = 1, 2, …… n, n is the number of frequencies, and n is a positive integer greater than or equal to 2; ，i The sparse linear system A = x b is obtained by discretizing the frequency-domain electromagnetic forward calculation through the finite element method or the finite difference method; The large-scale sparse linear system obtained by discretizing the same set of frequency-domain electromagnetic forward calculations through the finite element method (FEM) or the finite difference method (FDM): A x = b , where A is a sparse matrix, which is generated by the discretized partial differential equation and contains the physical parameters of the electromagnetic field (such as conductivity, permeability, etc.) and the discrete grid information; x is the vector to be solved, representing the values of the electromagnetic response characteristics to be solved of the exploration object (such as electric field or magnetic field) at the discrete grid nodes or elements; bis the right - hand vector, generated by source terms (such as current sources or magnetic field sources) and boundary conditions, representing external excitations or known conditions; different frequencies i of the sparse linear system is represented by A i x i = b i ( i = 1, 2, …… n), that is, the sparse linear system A x = b is transformed into the sparse linear system A of different frequencies i i x i = b i where A i is the sparse matrix of different frequencies, and A i is nrow × nrow complex square matrix; x i is the vector to be solved at different frequencies, with a dimension of nrow * 1 vector; b i is the right - hand - side vector at different frequencies, with a dimension of nrow * 1 vector ；i is the serial number at different frequencies; i = 1, 2, …… n, n is the number of frequencies, and n is a positive integer greater than or equal to 2; nrow is the number of degrees of freedom to be solved.

[0029] S11. Obtain the vector register length of the CPU processor; The vector lengths of current mainstream CPU processors are usually 256 or 512 bits, and most newly released CPUs in recent years are 512 bits.

[0030] Generally, the vector length will be stated on the CPU processor model, that is, the vector register length of the CPU processor, or by checking the CPU instruction set. For example, in the X86 architecture, if AVX256 is supported, it means the vector length >= 256; if AVX512 is supported, it means the vector length >= 512; this is prior art and will not be elaborated here.

[0031] S12. Calculate the number of vectors that a single machine instruction of the CPU can calculate simultaneously according to the vector register length and the calculation precision; Assume that the calculation precision adopted is 64 - bit floating - point and the vector register length of the computing platform is 512 bits. Then the number of vectors that a single machine instruction can calculate simultaneously is 512÷64 = 8.

[0032] S13. The sparse matrix A at different frequenciesi Convert to A i =(Re + Im i ), and then convert Im i to Im i = Im· f i ; where Re is the real matrix of the sparse matrix A i , and Im i is the imaginary matrix of the sparse matrix A i ; Im represents the reference complex matrix, f i is the frequency value; For different frequencies of the same set of meshes i the generated sparse matrix A i has the following characteristics: A i is a complex matrix and can be expressed as A i =(Re + Im i ), where Re and Im i represent the real matrix and the imaginary matrix respectively.

[0033] Matrices A with different frequencies i have exactly the same real part matrix Re; while the imaginary part matrix Im i has the following characteristics: Im i = Im· f i . Where Im represents the reference complex matrix, f i represents the frequency value, usually f i = 0.1, 2.0, 1.0, 10.0, 100.0, etc.

[0034] S14. Construct a frequency vector f i through the frequency value f vec , and f vec pad it to a multiple of the vector number. After padding, the vector length of the frequency vector f vec is vec_num; In step S12, it has been calculated that the number of vectorizations that can be calculated simultaneously by the current machine instruction is 8. Therefore, it is necessary to f vec pad it to a multiple of the vector number, and the value of the padded part can be any value in the frequency values. The length of the padded vector is vec_num.

[0035] For example, the number of frequencies to be solved is 5,f i = 0.1, 2.0, 1.0, 10.0, 100.0 need to be f i expanded into an array of length 8 and filled at the end, that is f i_vec = [0.1, 2.0, 1.0, 10.0, 100.0, f sup , f sup , f sup ]. f sup represents the filling value, f sup usually can take the maximum value of the frequency values participating in the calculation. In this embodiment f sup = 100.0, vec_num = 8. If the number of frequency numbers to be solved is 10, then vec_num = 16, and so on.

[0036] S15. Through the right - hand - side vector b i construct the right - hand - side b i_vec , b i_vec is nrow a vector of b i_vec rows, and each b i row is padded to a multiple of the number of vector elements, and the padded value is any element in the right - hand - side vector

[0037] Through the right - hand - side vector b i construct the right - hand - side b i_vec , b i_vec is nrow a vector of vec_num rows, and each element in it is a nrow - long vector,

[0038] Use b i [ j ] ( j = 0, 1, 2, 3 ……, nrow - 1), representing the i - th element of the right - hand - side corresponding to the frequency j After reconstructionb i_vec become b 0 j , b 1 j , b 2 j in a continuous manner. Similar to step S14, it is necessary to b i_vec perform vector padding, and the padding value can be selected from any b i array.

[0039] For example, if the number of frequencies to be solved is 5, then b i_vec the element arrangement order of n = nrow -1, b sup j represents the padding value.

[0040] {{ b 1 [0], b 2 [0], b 3 [0], b 4 [0], b 5 [0], b sup [0], b sup [0], b sup [0]}; { b 1 [1], b 2 [1], b 3 [1], b 4 [1], b 5 [1], b sup [1], b sup [1], b sup [1]}; … { b 1 [n],​​​​b 2 [n], b 3 [n], b 4 [n], b 5 [n], b sup [n], b sup [n], b sup [n]}}。

[0041] It can be understood that for the sparse linear system A i x i = b i After being converted into a vectorized calculation format, various variables can be stored for subsequent use in calculations. The specific storage method is not limited and can be stored in a cache, database, or file, which is selected according to specific needs.

[0042] S2. Initialization x i The approximate solution of x i_vec , and according to the maximum number of iterations max_iter and the convergence tolerance tolerance Use the vectorized calculation method to calculate the residual r i_vec and its two-norm β vec Then converge and update the approximate solution x i_vec ; S21. Set the maximum number of iterations max_iter and the convergence tolerance tolerance ; The maximum number of iterations max_iter and the convergence tolerance tolerance The values are determined according to actual needs; preferably, in this embodiment, the maximum number of iterations max_iter is set to 100, and the convergence tolerance tolerance is set to 10 -6 .

[0043] S22. Initialize the approximate solution x i_vec , x i_vec is a nrow row vector, and each element in it is a vec_num long vector, x i_vec The values in can usually take the value of 1; For example, if the number of frequencies to be solved is 5, then x i_vec the element arrangement order of nrow -1, x sup j represents a supplementary value.

[0044] {{ x 1 [0], x 2 [0], x 3 [0], x 4 [0], x 5 [0], x sup [0], x sup [0], x sup [0]}; { x 1 [1], x 2 [1], x 3 [1], x 4 [1], x 5 [1], x sup [1], x sup [1], x sup [1]}; …… { x 1 [n], x 2 [n], x 3 [n], x 4 [n], x 5 [n], x sup [n], x sup [n], x sup [n]}}.

[0045] S23. Vectorized calculation of the initial residual r ​i_vec , r i_vec is a nrow vector of length vec_num , and the calculation method is as follows; r i_vec = b i_vec -(Re + Im · f i_vec ) x i_vec ; In the process of calculating the product of a sparse matrix and a vector, the operation of multiplying each sparse matrix element by the vector is a multiplication of a vector by a vector; For example, when calculating (Re + Im · f i_vec ) x i_vec : (1) Take an element i, j of the matrix as an example; (2) Expand Re[i][j] into a vec_num vector Re vec [i][j]; (3) Calculate Im[i][j] · f i_vec into a vec_num vector Im vec [i][j]; (4) Obtain A vec [i][j] = Re vec [i][j] + Im vec [i][j]; (5) Calculate A vec [i][j] · x i_vec [j, :] = v vec , v vec is a vec_num vector.

[0046] The above calculations can all be executed through the vectorization instructions of modern CPUs. From the calculation method, it can be seen that each matrix multiplication can complete all matrix-vector multiplication operations related to frequency. The final r i_vec elements are arranged in the order shown below. Each row represents a vector, where n = nrow -1: {{ r 1 [0],r 2 [0], r 3 [0], r 4 [0], r 5 [0], r 6 [0], r 7 [0], r 8 [0]}; { r 1 [1], r 2 [1], r 3 [1], r 4 [1], r 5 [1], r 6 [1], r 7 [1], r 8 [1]}; …… { r 1 [n], r 2 [n], r 3 [n], r 4 [n], r 5 [n], r 6 [n], r 7 [n], r 8 [n]}};

[0047] S24. Calculate r i_vec of the two - norm β vec , β vec is a vector of dimension vec_num ; β vec = sqrt ( r 0_vec 2 + r 1_vec2 + r 2_vec 2 +…… r n_vec 2 ); wherein r 0_vec represents a set of vectors, { r 1 [0], r 2 [0], r 3 [0], r 4 [0], r 5 [0], r 6 [0], r 7 [0], r 8 [0]}.

[0048] The above calculation can be performed through the vectorized calculation instructions of a modern CPU.

[0049] Corresponding to the calculated β vec [i]= sqrt ( r i [0] 2 + r i [1] 2 + r i [2] 2 +…… r i [n] 2 ) (i = 1, 2... 8); S25. Convergence judgment. If β vec each value in tolerance satisfies being less than the convergence tolerance x i_vec .

[0050] S3. Initialize the orthogonal basis matrix V vec and the Hessenberg matrix H vec , and calculate V vec [:, 0], V vec [:,0] = ri_vec / β vec ; S31. Initialize the orthogonal basis matrix V vec , where the orthogonal basis matrix V vec is a matrix of dimension nrow ×( max_iter + 1), and each matrix element is a vector of length vec_num for storing the orthogonal basis of the Krylov subspace. Initialize the value of V vec to 0; S32. Calculate V vec [:, 0], V vec [:, 0] = r i_vec / β vec ; V vec [:, 0] stores the results as follows. Each row represents a vector, where n = nrow -1: {{ r 1 [0] / β [0], r 2 [0] / β [1], …… r 8 [0] / β [7]}; { r 1 [1] / β [0], r 2 [1] / β [1], …… r 8 [1] / β [7]}; …… { r 1 [n] / β [0], r 2 [n] / β [1], …… r 8 [n] / β [7]}};

[0051] In GMRES (Generalized Minimal Residual Method), the orthogonal basis matrix V vec is a matrix composed of a set of orthogonal basis vectors generated by the Arnoldi process. The orthogonal basis matrix V vec is an auxiliary matrix constructed by GMRES to solve the sparse linear system.

[0052] S33. Initialize the Hessenberg matrix H vec , the matrix H vec constructs a matrix with dimensions of max_iter × ( max_iter +1), and each matrix element is a vector of length vec_num . The value of H vec is initialized to 0.

[0053] The said Hessenberg matrix H vec is an upper Hessenberg matrix; S4. Update the orthogonal basis matrix V vec and the Hessenberg matrix H vec through iterative orthogonalization, Hessenberg matrix update, and least squares solution, and calculate the current residual Res vec . When the current residual Res vec converges or the number of iterations reaches the maximum, calculate the updated approximate solution V vec according to the orthogonal basis matrix x i_vec , so as to obtain the solution of the frequency-domain electromagnetic forward calculation.

[0054] In step S4, it further includes: obtaining the electromagnetic response characteristics of the exploration object according to the solution of the frequency-domain electromagnetic forward calculation.

[0055] This step mainly performs the main loop iteration: S41. Set the current iteration loop count iter = 0; S42. Calculate the orthogonal vector w i_vec = (Re + Im· f i_vec ) * V vec [:, iter , w i_vec is a vector with nrow rows, and each element in it is a vector of length vec_num ; The matrix-vector multiplication operation involved can refer to the process of step S23 in step S2; S43. Modified orthogonalization, the specific steps are as follows: S431. Set k = 0; Enter the loop; S432. Projection coefficient H vec [k] iter Calculate: H vec [k] iter = dot (( V vec [:, k ) T , w i_vec ), dot indicating the sum of the dot products of two vectors; S433. Through w i_vec = w i_vec -H vec k iter V vec [:, k for orthogonalization correction; S434. If k is less than iter + 1 ,k = k + 1, jump to S432, otherwise jump to S435; S435. Modified orthogonalization ends; The matrix-vector multiplication operation involved can refer to the process of step S23 in step S2; S44. Calculate the 2-norm of the orthogonal vector w i_vec and assign it to the Hessenberg matrix H vec ; H vec iter + 1] iter = sqrt ( w 0_vec 2 + w 1_vec 2 + w 2_vec 2 + …… w n_vec 2 ), where n = nrow - 1; ​​​​​​For the detailed calculation method, refer to the calculation of the vector two-norm in step S24 of step S2.

[0056] S45. Normalize the Hessenberg matrix H vec is an orthogonal basis matrix V vec ; V vec [:, iter+1] = w i_vec / H vec iter +1] iter ; The process between S42 and S45 is the process of orthogonal basis expansion; S46. Solve the Hessenberg matrix H vec by the least squares method to obtain the first vector array y i_vec and the second vector array e i_vec ; This part cannot be vectorized and can only be solved sequentially; the solution process is as follows: S461. Set k =0; Enter the loop; S462. Construct H k as a scalar matrix of ( iter +2) × ( iter +1), and H k stores the vec [: iter +2, : iter +1] in each vector element of the k +1-th element; For example, H 0 stores the vec [: iter +2, : iter +1] in each vector element of the 1st element; S463. Construct e k as a scalar array of ( iter +2), initialize e k to 0; e k [0] = β k ; S464. Solve H k y k =e k , and obtain y ​​k is a scalar array of ( iter + 1); S465. If k is less than vec_num,k = k +1 ， jump to S462; otherwise, jump to S466; S466. Assemble y k (k = 0, 1,..., vec_num - 1) into a new first vector array y i_vec , y i_vec which is a iter + 1-row vector, and each element of which is a vec_num -long vector; S467. Assemble e k (k = 0, 1,..., vec_num - 1) into a new second vector array e i_vec , e i_vec which is a iter + 2-row vector, and each element of which is a vec_num -long vector.

[0057] S47. Calculate the current residual Res vec based on the Hessenberg matrix H y i_vec , the first vector array e i_vec and the second vector array vec and check for convergence; S471. Calculate tmp vec = H vec [: iter + 2, : iter + 1] y i_vec -e i_vec , tmp vec which is a iter + 2-row vector, and each element of which is a vec_num -long vector.

[0058] S472. Calculate Res vec for the current iteration. The residual value of the current iteration is as follows. Res vec is a vector of length vec_num: Resvec = sqrt ( tmp 0_vec 2 + tmp 1_vec 2 + tmp 2_vec 2 +…… tmp n_vec 2 ); S473. Determine whether all values in the current residual Res vec are less than the convergence tolerance tolerance . If so, jump to step S49; if not, jump to step S48; S48. Determine whether iter is less than the maximum number of iterations max_iter. If so, jump to step S42; if not, jump to S49; S49. Calculate the final solution and output the result.

[0059] x i_vec = x i_vec +V vec [:, iter + 1] y i_vec。

[0060] The obtained x i_vec solution in is the solution of the frequency-domain electromagnetic forward modeling. The storage information has been given in step S22 of S2, as follows. Then x i_vec The x 1 [0], x 1 [1]…… x 1 [n] corresponds to the frequency of f 1 Similarly x i [0], x i [1]…… x i [n] corresponds to the frequency of f i The supplementary part of the value can be discarded.

[0061] {{ x 1 [0], x 2 [0],x 3 [0], x 4 [0], x 5 [0], x sup [0], x sup [0], x sup [0]}; { x 1 [1], x 2 [1], x 3 [1], x 4 [1], x 5 [1], x sup [1], x sup [1], x sup [1]}; …… { x 1 [n], x 2 [n], x 3 [n], x 4 [n], x 5 [n], x sup [n], x sup [n], x sup [n]}}.

[0062] Finally, obtain the electromagnetic response characteristics of the exploration object according to the solution of the frequency-domain electromagnetic forward calculation.

[0063] In the prior art, magnetotelluric method (MT) and controlled-source electromagnetic method (CSEM) both belong to the frequency-domain electromagnetic method and are widely used in fields such as metal ore exploration, oil and gas exploration, and geothermal resource detection. Frequency-domain electromagnetic forward calculation is one of the core technologies of the frequency-domain electromagnetic method. Traditional solution methods need to solve for frequencies sequentially, with low speed and efficiency; or solve through multi-core parallel computing, which requires additional computing resources.

[0064] During metal ore exploration, metal ore bodies usually exhibit significant electrical anomalies (such as low resistivity or high resistivity), but their burial depth, shape, and scale are complex and variable. Traditional forward modeling methods are time-consuming in calculation and difficult to meet the needs of large-scale exploration. The vectorized parallel solution method for frequency-domain electromagnetic forward calculation in this embodiment can greatly improve the data processing speed, quickly identify electromagnetic response characteristics such as the position of the ore body and electrical parameters, shorten the data processing cycle by more than 50%, significantly reduce the hardware cost, improve the accuracy of ore body boundary identification at the same time, reduce the blindness of drilling, and lower the exploration cost.

[0065] During oil and gas exploration, oil and gas reservoirs usually show high resistivity anomalies, but their electromagnetic responses are affected by complex formation structures (such as salt domes, faults). Traditional forward modeling methods are difficult to quickly simulate electromagnetic response characteristics such as the electromagnetic field distribution at multiple frequency points and multiple scales. The vectorized parallel solution method for frequency-domain electromagnetic forward calculation in this embodiment can reduce the dependence on high-performance computing clusters, efficiently simulate electromagnetic responses under complex geological conditions, shorten the data processing cycle by more than 50%, significantly reduce the hardware cost, improve the accuracy of oil and gas reservoir boundary identification at the same time, and reduce the exploration risk.

[0066] This embodiment can make full use of the long vector register characteristics of modern CPU processors, can solve multiple frequency-related sparse linear systems simultaneously at one time, can greatly improve the calculation speed, and can significantly improve the efficiency of frequency-domain electromagnetic forward calculation.

[0067] This embodiment proposes a vectorized parallel solution method for frequency-domain electromagnetic forward calculation, including: S1. Convert sparse linear systems A at different frequencies into vectorized calculation formats according to the vector register length and calculation accuracy of the CPU processor; i x i = b i S2. Initialize the approximate solution of, and calculate the residual and its two-norm using the vectorized calculation method, and then converge and update the approximate solution; x i of x i_vec S3. Initialize the orthogonal basis matrix and the Hessenberg matrix H, and calculate [:, 0], [:, 0] = / r i_vec and its two-norm β vec and then converge and update the approximate solution x i_vec ; V vec and the Hessenberg matrix H vec and calculate V vec [:, 0], V vec [:, 0] = r i_vec / β vec ; S4. Update the orthogonal basis matrix V vec and the Hessenberg matrix H vec through iterative orthogonalization, Hessenberg matrix update, and least squares solution, and calculate the current residual Res vec , and when the current residual Res vec converges or the maximum number of iterations is reached, calculate the updated approximate solution V vec according to the orthogonal basis matrix x i_vec . In this embodiment, the left-end matrix and the right-end term of the large-scale sparse linear system obtained after discretization in the frequency-domain electromagnetic forward calculation process are converted into a vectorized storage format. Through the idea of batch processing, sparse linear systems at multiple frequencies are solved simultaneously, and calculations are performed through the vector instructions of modern CPUs. When using the generalized minimum residual method to solve sparse linear systems, each core operator in it (mainly including calculations such as sparse matrix-vector multiplication, vector dot multiplication, and vector scalar multiplication) is extended into a frequency-dependent vector calculation, and a convergence criterion based on vectorized calculation is set to achieve the goal of solving sparse linear systems at multiple frequencies at one time; this embodiment combines the calculation characteristics and memory access characteristics of multi-frequency sparse linear systems, makes full use of the long vector registers of the CPU, reduces the memory access times of the matrix, solves the problems of reduced calculation efficiency and low utilization rate of calculation resources in the prior art, and realizes fast and efficient multi-frequency frequency-domain electromagnetic forward calculation, greatly improving the calculation efficiency of frequency-domain electromagnetic forward calculation.

[0068] In addition, compared with traditional frequency parallel computing, this embodiment does not require any additional computing resources, and the saved CPU cores can be used for other calculations, such as parallel computing of matrix rows, or calculations of other programs, which can save more computing resources.

[0069] Embodiment 2 Referring to Figure 2 , this embodiment provides a vectorized parallel solution device for frequency-domain electromagnetic forward calculation, which includes the following units: A vectorization unit, configured to convert sparse linear systems A at different frequencies i x i = b i into a vectorized calculation format according to the vector register length and calculation accuracy of the CPU processor, including the real matrix Re of the sparse matrix A i , the imaginary matrix Im of the sparse matrix A i , the reference complex matrix Im, and the frequency value i . fi , frequency vectors f vec and the right - hand side terms b i_vec ; where, A i is a sparse matrix of different frequencies, and A i is nrow × nrow complex square matrix; x i is the vector to be solved for different frequencies, and the dimension is nrow * 1 vector; b i is the right - hand side vector for different frequencies, and the dimension is nrow * 1 vector ；i is the sequence number for different frequencies; i = 1, 2, …… n, n is the number of frequencies, and n is a positive integer greater than or equal to 2; nrow is the number of degrees of freedom to be solved; The first initialization unit is used to initialize x i approximate solution of x i_vec , and calculate the residual r i_vec and its two - norm β vec and then converge to update the approximate solution x i_vec ; The second initialization unit is used to initialize the orthogonal basis matrix V vec and the Hessenberg matrix H vec , and calculate V vec [:, 0], V vec [:, 0] = r i_vec / β vec ; The iterative calculation unit is used to update the orthogonal basis matrix V vec and the Hessenberg matrix H vec through iterative orthogonalization, Hessenberg matrix update, and least - squares solution, and calculate the current residual Res vec , and when the current residual Res vec converges or the maximum number of iterations is reached, calculate and update the approximate solution V vec according to the orthogonal basis matrix x i_vec, thus obtaining the solution of the frequency-domain electromagnetic forward calculation.

[0070] Embodiment 3 Reference Figure 3 , Figure 3 is a schematic structural diagram of a vectorized parallel solution device for frequency-domain electromagnetic forward calculation in this embodiment. The vectorized parallel solution device 20 for frequency-domain electromagnetic forward calculation in this embodiment includes a processor 21, a memory 22, and a computer program stored in the memory 22 and executable on the processor 21. When the processor 21 executes the computer program, the steps in the above method embodiment are implemented. Alternatively, when the processor 21 executes the computer program, the functions of each module / unit in the above device embodiments are implemented.

[0071] Exemplarily, the computer program can be divided into one or more modules / units. The one or more modules / units are stored in the memory 22 and executed by the processor 21 to complete the present invention. The one or more modules / units can be a series of computer program instruction segments capable of performing specific functions, and these instruction segments are used to describe the execution process of the computer program in the vectorized parallel solution device 20 for frequency-domain electromagnetic forward calculation. For example, the computer program can be divided into the respective modules in Embodiment 2. For the specific functions of each module, please refer to the working process of the device described in the above embodiments, and details are not described herein again.

[0072] The vectorized parallel solution device 20 for frequency-domain electromagnetic forward calculation may include, but is not limited to, a processor 21 and a memory 22. Those skilled in the art can understand that the schematic diagram is only an example of the vectorized parallel solution device 20 for frequency-domain electromagnetic forward calculation, and does not constitute a limitation on the vectorized parallel solution device 20 for frequency-domain electromagnetic forward calculation. It may include more or fewer components than shown in the figure, or combine certain components, or different components. For example, the vectorized parallel solution device 20 for frequency-domain electromagnetic forward calculation may further include an input / output device, a network access device, a bus, etc.

[0073] The processor 21 may be a Central Processing Unit (CPU), or may also be other general-purpose processors, Digital Signal Processors (DSPs), Application Specific Integrated Circuits (ASICs), Field-Programmable Gate Arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or the processor may also be any conventional processor, etc. The processor 21 is the control center of the vectorized parallel solving device 20 for frequency-domain electromagnetic forward calculation, and connects various parts of the vectorized parallel solving device 20 for frequency-domain electromagnetic forward calculation through various interfaces and lines.

[0074] The memory 22 can be used to store the computer programs and / or modules. The processor 21 realizes various functions of the vectorized parallel solving device 20 for frequency-domain electromagnetic forward calculation by running or executing the computer programs and / or modules stored in the memory 22, and by calling the data stored in the memory 22. The memory 22 mainly includes a program storage area and a data storage area. Among them, the program storage area can store an operating system, application programs required for at least one function (such as a sound playback function, an image playback function, etc.); the data storage area can store data created according to the use of the mobile phone (such as audio data, phone book, etc.). In addition, the memory 22 may include high-speed random access memory, and may also include non-volatile memory, such as a hard disk, a memory, a plug-in hard disk, a Smart Media Card (SMC), a Secure Digital (SD) card, a Flash Card, at least one magnetic disk storage device, a flash memory device, or other volatile solid-state storage devices.

[0075] Among them, if the modules / cells integrated in the vectorized parallel solving device 20 for frequency-domain electromagnetic forward calculation are implemented in the form of software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on such understanding, to implement all or part of the processes in the above-described embodiment methods of the present invention, it can also be completed by a computer program instructing relevant hardware. The computer program can be stored in a computer-readable storage medium. When the computer program is executed by the processor 21, the steps of the above-described method embodiments can be implemented. Among them, the computer program includes computer program code, and the computer program code can be in the form of source code, object code, executable file, or some intermediate form, etc. The computer-readable medium can include: any entity or device capable of carrying the computer program code, recording medium, USB flash drive, mobile hard disk, magnetic disk, optical disc, computer memory, read-only memory (ROM), random access memory (RAM), electrical carrier signal, telecommunication signal, and software distribution medium, etc. It should be noted that the content included in the computer-readable medium can be appropriately increased or decreased according to the requirements of legislation and patent practice in the jurisdiction. For example, in some jurisdictions, according to legislation and patent practice, the computer-readable medium does not include electrical carrier signals and telecommunication signals.

[0076] It should be noted that the device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separated, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed to multiple network units. Some or all of the modules can be selected according to actual needs to achieve the purpose of the solution of this embodiment. In addition, in the attached drawings of the device embodiments provided by the present invention, the connection relationship between the modules indicates that they have a communication connection, which can be specifically implemented as one or more communication buses or signal lines. Those of ordinary skill in the art can understand and implement it without creative effort.

[0077] The parts not detailed in the present invention are prior art. For those skilled in the art, it is obvious that the present invention is not limited to the details of the above-described exemplary embodiments, and can be implemented in other specific forms without departing from the spirit or basic characteristics of the present invention. Therefore, from any point of view, the embodiments should be regarded as exemplary and non-limiting, aiming to include all changes falling within the meaning and scope of the equivalent elements in the present invention.

Claims

1. A vectorized parallel solution method for frequency domain electromagnetic forward modeling, characterized in that: The specific steps include: S1. According to the vector register length and calculation accuracy of the CPU processor, the sparse linear system A with different frequencies is i x i = b i Convert to vectorized computing format, including sparse matrix A i The real matrix Re and the sparse matrix A i The imaginary matrix Im i , reference complex matrix Im, frequency value f i , frequency vector f vec and right end b i_vec ; Among them, A i is a sparse matrix of different frequencies, A i for nrow × nrow The complex matrix of ; x i is the vector to be solved at different frequencies, with dimensions of nrow*1 vector of; b i is the right-hand term vector of different frequencies, with dimension nrow*1 Vector ；i are the serial numbers of different frequencies; i =1,2,……n, n is the number of frequencies, n is a positive integer greater than or equal to 2; nrow is the number of degrees of freedom that need to be solved; S2. Initialization x i The approximate solution of x i_vec , and use vectorized calculation method to calculate the residual r i_vec and its bi-norm β vec Post-convergence update of approximate solution x i_vec ; S3. Initialize the orthogonal basis matrix V vec and the Hessenberg matrix H vec , and calculate V vec [:, 0], V vec [:, 0] = r i_vec / β vec ; S4. Orthogonal basis matrix by iterative orthogonalization, Hessenberg matrix update and least squares solution V vec and the Hessenberg matrix H vec Update and calculate the current residual Res vec , and in the current residual Res vec When convergence or the number of iterations is maximum, the orthogonal basis matrix V vec Compute updated approximate solution x i_vec , thus obtaining the solution of the frequency domain electromagnetic forward modeling calculation.

2. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 1 is characterized in that: Step S1 specifically includes the following steps: S11, obtaining the vector register length of the CPU processor; S12, calculating the number of vectors that can be calculated simultaneously by one machine instruction of the CPU according to the vector register length and calculation accuracy; S13, the sparse matrix A of different frequencies i Convert to A i =(Re+Im i ), then Im i Convert to Im i = Im f i ; Where Re is the sparse matrix A i The real matrix, Im i A is a sparse matrix i Im represents the base complex matrix. f i is the frequency value; S14, passing frequency value f i Constructing the frequency vector f vec , and f vec Pad to a multiple of the number of vectors, the frequency vector after padding f vec The vector length is vec_num; S15, through the right end vector b i Construct the right hand side b i_vec , b i_vec for nrow The vector of rows and b i_vec Each row is padded to a multiple of the number of vectors, and the padded value is the right end vector b i Any element in .

3. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 2 is characterized in that: Step S2 includes the following steps: S21. Set the maximum number of iterations max_iter and convergence tolerance tolerance ; S22, Initialization x i The approximate solution of x i_vec , x i_vec For one nrow A vector of rows, each element of which is a vec_num Long vector; S23, vectorized calculation of initial residual r i_vec : r i_vec = b i_vec -(Re+ Im· f i_vec ) x i_vec ; r i_vec For one nrow A vector of rows, each element of which is a vec_num Long vector; S24. Calculate initial residual r i_vec The second norm of β vec : β vec = sqrt ( r 0_vec 2 + r 1_vec 2 + r 2_vec 2 +…… r n_vec 2 )； in r n_vec represents a set of vectors, β vec For a dimension vec_num Vector, n= nrow -1.

4. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 3 is characterized in that: Step S3 specifically includes the following steps: S31, initialize the orthogonal basis matrix V vec , the orthogonal basis matrix V vec For Dimension nrow ×( max_iter +1) matrix, each matrix element is a vec_num vector of lengths; S32, calculation V vec [:, 0] ，V vec [:, 0] = r i_vec / β vec ; S33. Initialize the Hessenberg matrix H vec , the matrix H vec The construction dimension is max_iter ×( max_iter +1) matrix, each matrix element is a vec_num Vector of lengths.

5. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 4 is characterized in that: Step S4 specifically includes the following steps: S41. Set the current number of iteration cycles iter =0; S42. Calculate orthogonal vectors w i_vec = (Re + Im f i_vec )* V vec [: , iter ], w i_vec is a nrow A vector of rows, each element of which is a vec_num Long vector; S43. Modify orthogonalization; S44. Calculate orthogonal vectors w i_vec The second norm of and assign it to the Hessenberg matrix H vec : H vec [ iter +1][ iter ] =sqrt ( w 0_vec 2 + w 1_vec 2 + w 2_vec 2 +…… w n_vec 2 ), where n= nrow -1; S45. Standardized Hessenberg matrix H vec is an orthogonal basis matrix V vec : V vec [:, iter+1] = w i_vec / H vec [ iter +1][ iter ]; S46. Use the least squares method to calculate the Hessenberg matrix H vec Solve to get the first vector array y i_vec and the second vector array e i_vec ; S47, according to the Hessenberg matrix H vec , first vector array y i_vec and the second vector array e i_vec Calculate the current residual Res vec And check whether it converges; S48, if the current residual Res vec If it does not converge, determine whether iter is less than the maximum number of iterations max_iter. If so, jump to step S42; if not, jump to step S49; S49, according to the orthogonal basis matrix V vec and the first vector array y i_vec For approximate solution x i_vec Update to obtain the solution of frequency domain electromagnetic forward modeling: x i_vec = x i_vec +V vec [:,iter+1] y i_vec .

6. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 5 is characterized in that: Step S43 specifically includes the following steps: S431, Settings k =0; Enter the loop; S432, projection coefficient H vec [k][ iter ]calculate: H vec [k][ iter ] =dot (( V vec [: , k ] ) T , w i_vec ), dot Represents the sum of two vector dot products; S433, through w i_vec = w i_vec - H vec [ k ][ iter ] V vec [: , k ] to perform orthogonal correction; S434, if k Less than iter+1 , k=k+1, Jump to S432, otherwise jump to S435; S435, the corrected orthogonalization is completed.

7. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 5 is characterized in that: Step S46 specifically includes the following steps: S461, Settings k =0; Enter the loop; S462, Structure H k for( iter +2)×( iter +1), H k Storage H vec [: iter +2, : iter +1] in each vector element k +1 element; S463, structure e k for( iter +2) scalar array, initialize e k is 0; e k [0] = β [ k ]; S464, solve H by least squares method k y k =e k ,get y k is a ( iter +1) scalar array; S465, if k Less than vec_num, k = k +1 ， Jump to S462, otherwise jump to S466; S466, will y k (k=0, 1, ... vec_num -1) Assemble into a new first vector array y i_vec , y i_vec is a iter +1 row vector, each element of which is a vec_num Long vector; S467, will e k (k=0, 1, ... vec_num -1) Assemble into a new second vector array e i_vec , e i_vec is a iter +2 rows of vectors, each element of which is a vec_num Long vector.

8. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 5, characterized in that: Step S47 specifically includes the following steps: S471, Calculation tmp vec =H vec [: iter +2, : iter +1] y i_vec -e i_vec , tmp vec is a iter +2 rows of vectors, each element of which is a vec_num Long vector; S472, Calculate Res vec , the residual value of the current iteration is as follows, Res vec is a vector of length vec_num: Res vec = sqrt ( tmp 0_vec 2 + tmp 1_vec 2 + tmp 2_vec 2 +…… tmp n_vec 2 ); S473, determine the current residual Res vec Are all values ​​in less than the convergence tolerance tolerance If yes, jump to step S49, if no, jump to step S48.

9. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 1, characterized in that: The sparse linear system A x = b It is obtained by discretizing the frequency domain electromagnetic forward modeling through the finite element method or the finite difference method.

10. A vectorized parallel solution device for frequency domain electromagnetic forward modeling, characterized in that: The following units are included: The vectorization unit is used to convert the sparse linear system A with different frequencies into i x i = b i Convert to vectorized computing format, including sparse matrix A i The real matrix Re and the sparse matrix A i The imaginary matrix Im i , reference complex matrix Im, frequency value f i , frequency vector f vec and right end b i_vec ; Among them, A i is a sparse matrix of different frequencies, A i for nrow × nrow The complex matrix of ; x i is the vector to be solved at different frequencies, with dimensions of nrow*1 vector of; b i is the right-hand term vector of different frequencies, with dimension nrow*1 Vector ；i are the serial numbers of different frequencies; i =1,2,……n, n is the number of frequencies, n is a positive integer greater than or equal to 2; nrow is the number of degrees of freedom that need to be solved; The first initialization unit is used to initialize x i The approximate solution of x i_vec , and use vectorized calculation method to calculate the residual r i_vec and its bi-norm β vec Post-convergence update of approximate solution x i_vec ; The second initialization unit is used to initialize the orthogonal basis matrix V vec and the Hessenberg matrix H vec , and calculate V vec [:,0], V vec [:, 0] = r i_vec / β vec ; Iterative calculation unit, used to solve the orthogonal basis matrix through iterative orthogonalization, Hessenberg matrix update and least squares solution V vec and the Hessenberg matrix H vec Update and calculate the current residual Res vec , and in the current residual Res vec When convergence or the number of iterations is maximum, the orthogonal basis matrix V vec Compute updated approximate solution x i_vec .

Images