A vectorized parallel solution method and device for frequency domain electromagnetic forward calculation

Abstract

The present invention provides a vectorized parallel solution method and device for frequency domain electromagnetic forward modeling, which relates to the field of algorithm solution technology. The present invention converts the large-scale sparse linear system obtained after discretization in the frequency domain electromagnetic forward modeling process into a vectorized storage format, and solves the sparse linear systems of multiple frequencies at the same time by using the vector instructions of modern CPUs and the idea of ​​batch processing. When using the generalized minimum residual method to solve the sparse linear system, each core operator is expanded into a frequency-related vector calculation to achieve the goal of solving the sparse linear systems of multiple frequencies at one time; the present invention combines the computational characteristics and memory access characteristics of multi-frequency sparse linear systems, fully utilizes the long vector registers of the CPU, reduces the number of memory accesses of the matrix, solves the problems of reduced computational efficiency and low computational resource utilization in the existing technology, realizes fast and efficient multi-frequency frequency domain electromagnetic forward modeling, and greatly improves computational efficiency.

Description

Technical Field

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

[0002] Frequency-domain electromagnetic forward modeling is one of the core technologies in geophysical exploration (especially electromagnetic methods). Frequency-domain electromagnetic forward modeling can simulate the propagation behavior of electromagnetic waves in underground media at different frequencies of exploration objects, thereby obtaining key information about the resistivity distribution, electrical layer thickness, geometry, electrical interface position, phase information, and the size and location of anomalies of the underground structure, thereby obtaining the electromagnetic response characteristics of the exploration object. These electromagnetic response characteristics are of great significance to fields such as geological exploration, resource surveys, and environmental monitoring.

[0003] Conventional forward modeling, after discretization using the finite element method (FEM) or finite-difference time-domain (FDTD) method, typically generates large, sparse linear systems. These systems typically account for the majority of the forward modeling computation time, often exceeding 80%. The generalized minimum residual method (GMRES) is an iterative algorithm for efficiently solving large, sparse, asymmetric linear systems and is widely used in electromagnetic field forward modeling. Computational performance and computational time are crucial to actual production efficiency. Therefore, it is necessary to improve forward modeling computation speed by leveraging the architectural characteristics of modern general-purpose CPUs.

[0004] The commonly used calculation methods for frequency domain electromagnetic forward modeling are: one is to solve the linear system generated by frequencies sequentially; the other is to solve in parallel by distributing different frequencies on multiple CPU cores.

[0005] Current CPU processors typically provide 256-512-bit vector registers, capable of simultaneously performing 4-8 64-bit floating-point calculations or 8-16 32-bit floating-point calculations. However, commonly used computational methods fail to incorporate the computational and memory access characteristics of multi-frequency sparse linear systems. This leads to repeated reading of sparse linear matrices and inability to fully utilize the long vectorization capabilities of current CPU processors, resulting in reduced computational efficiency and low resource utilization.

[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 domain electromagnetic forward modeling calculations are difficult problems that need to be solved urgently. Summary of the Invention

[0007] In order to overcome the deficiencies in the background technology, the present invention provides a vectorized parallel solution method and apparatus for frequency domain electromagnetic forward modeling.

[0008] In order to achieve the above-mentioned object of the invention, the present invention adopts the following technical solutions:

[0009] In a first aspect, the present invention provides a vectorized parallel solution method for frequency domain electromagnetic forward modeling, comprising the following steps:

[0010] S1, according to the vector register length and calculation accuracy of the CPU processor, the sparse linear system A of different frequencies is converted into i x i = b i Convert to vectorized calculation format, including sparse matrix A i The real matrix Re and sparse matrix A i The imaginary matrix Im i , reference complex matrix Im, frequency value f i , frequency vector f vec and the 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; b i is the right-hand term vector of different frequencies, with dimension nrow * 1 Vector ；i is the serial number 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;

[0011] S2. Initialization x i Approximate solution of x i_vec , and use vectorized calculation method to calculate the residual r i_vec and its second norm β vec Post-convergence update of approximate solution x i_vec ;

[0012] 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 ;

[0013] S4, 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 Calculate updated approximate solution x i_vec .

[0014] Specifically, step S1 includes the following steps:

[0015] S11. Obtain the vector register length of the CPU processor;

[0016] S12. Calculate the number of vectors that can be simultaneously calculated by one CPU machine instruction based on the vector register length and calculation accuracy;

[0017] 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 is a sparse matrix A i Im represents the base complex matrix, f i is the frequency value;

[0018] S14, passing frequency value f i Constructing 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;

[0019] S15, through the right end vector bi 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 .

[0020] Specifically, step S2 includes the following steps:

[0021] S21. Set the maximum number of iterations max_iter and convergence tolerance tolerance ;

[0022] S22, Initialization x i 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;

[0023] S23, vectorized calculation of initial residual r i_vec :

[0024] r i_vec = b i_vec -(Re+ Im· f vec ) x i_vec ;

[0025] r i_vec For one nrow A vector of rows, each element of which is a vec_num long vector;

[0026] S24. Calculate initial residual r i_vec The second norm of β vec :

[0027] β vec = sqrt ( r 0_vec 2 + r 1_vec 2 +r 2_vec 2 +…… r n_vec 2 );

[0028] in r i_vec represents a set of vectors, β vec For a dimension vec_num vector, n= nrow -1.

[0029] Specifically, step S3 includes the following steps:

[0030] S31. Initialize the orthogonal basis matrix V vec , the orthogonal basis matrix V vec Dimension nrow ×( max_iter +1) matrix, each matrix element is a vec_num vector of lengths;

[0031] S32, calculation V vec [:, 0], V vec [:, 0] = r i_vec / β vec ;

[0032] S33. Initialize the Hessenberg matrix H vec , matrix H vec The construction dimension is max_iter ×( max_iter +1) matrix, each matrix element is a vec_num A vector of lengths.

[0033] Specifically, step S4 includes the following steps:

[0034] S41. Set the current number of iteration cycles iter =0;

[0035] S42. Calculate orthogonal vectors w i_vec = (Re + Im· f vec ) * V vec [: , iter ], w i_vec is anrow A vector of rows, each element of which is a vec_num long vector;

[0036] S43. Modify orthogonalization;

[0037] S44. Calculate orthogonal vectors w i_vec The second norm of and assign it to the Hessenberg matrix H vec :

[0038] 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;

[0039] S45. Normalized Hessenberg matrix H vec is an orthogonal basis matrix V vec :

[0040] V vec [:, iter+1] = w i_vec / H vec [ iter +1][ iter ];

[0041] S46, through the least squares method to Hessenberg matrix H vec Solve to get the first vector array y i_vec and the second vector array e i_vec ;

[0042] 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;

[0043] 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; otherwise, jump to step S49.

[0044] S49, according to the orthogonal basis matrix V vec and the first vector array y i_vec 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 .

[0045] Specifically, step S43 includes the following steps:

[0046] S431, Settings k =0; Enter the loop;

[0047] S432, projection coefficient H vec [k][ iter ]calculate:

[0048] H vec [k][ iter ] = dot (( V vec [: , k ] ) T , w i_vec ), dot Represents the sum of two vector dot products;

[0049] S433, through w i_vec = w i_vec -H vec [ k ][ iter ] V vec [: , k ] to perform orthogonal correction;

[0050] S434, if k Less than iter+1 ,k = k + 1, Jump to S432, otherwise jump to S435;

[0051] S435. Correction orthogonalization ends.

[0052] Specifically, step S46 includes the following steps:

[0053] S461, Settings k =0; Enter the loop;

[0054] S462, Structure H k for( iter +2)×( iter +1) scalar matrix, H k Storage H vec [: iter +2, : iter +1] in each vector element k +1 element;

[0055] S463, structure e k for( iter +2) scalar array, initialize e k is 0; e k [0] = β [ k ];

[0056] S464, solve H by least squares method k y k =e k ,get y k is a ( iter +1);

[0057] S465, if k Less than vec_num,k = k +1 ， Jump to S462, otherwise jump to S466;

[0058] 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;

[0059] S467, will e k (k=0, 1, ... vec_num -1) Assemble into a new second vector array ei_vec , e i_vec is a iter +2 rows of vectors, each element of which is a vec_num Long vector.

[0060] Specifically, step S47 includes the following steps:

[0061] S471, calculation tmp vec =H vec [: iter +2, : iter +1] y i_vec -e i_vec ,

[0062] tmp vec is a iter +2 rows of vectors, each element of which is a vec_num long vector;

[0063] S472, calculate Res vec , the residual value of the current iteration is as follows, Res vec is a vector of length vec_num:

[0064] Res vec = sqrt ( tmp 0_vec 2 + tmp 1_vec 2 + tmp 2_vec 2 +…… tmp n_vec 2 );

[0065] S473, judging the current residual Res vec Are all values ​​in less than the convergence tolerance tolerance If yes, go to step S49, if no, go to step S48.

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

[0067] Specifically, step S4 further includes: acquiring the electromagnetic response characteristics of the exploration object according to the solution of the frequency domain electromagnetic forward calculation.

[0068] In a second aspect, the present invention provides a vectorized parallel solution device for frequency domain electromagnetic forward modeling, comprising the following units:

[0069] The vectorization unit is used to convert the sparse linear system A of different frequencies into i x i = b i Convert to vectorized calculation format, including sparse matrix A i The real matrix Re and sparse matrix A i The imaginary matrix Im i , reference complex matrix Im, frequency value f i , frequency vector f vec and the 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; b i is the right-hand term vector of different frequencies, with dimension nrow * 1 Vector ；i is the serial number 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;

[0070] The first initialization unit is used to initialize x i Approximate solution of x i_vec , and use vectorized calculation method to calculate the residual r i_vec and its second norm β vec Post-convergence update of approximate solution x i_vec ;

[0071] 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], Vvec [:, 0] = r i_vec / β vec ;

[0072] 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 Calculate updated approximate solution x i_vec , thus obtaining the solution of the frequency domain electromagnetic forward calculation.

[0073] In a third aspect, the present invention also provides an electronic device comprising a processor, a memory, a communication interface, and one or more programs, wherein the one or more programs are stored in the memory and configured to be executed by the processor, and the program includes instructions for executing the steps in the method described in any one of the first aspects.

[0074] The present invention proposes a vectorized parallel solution method for frequency domain electromagnetic forward modeling, comprising: S1, according to the vector register length and calculation accuracy of the CPU processor, the sparse linear system A of different frequencies is converted into i x i = b i Convert to vectorized calculation format; S2, initialization x i Approximate solution of x i_vec , and use vectorized calculation method to calculate the residual r i_vec and its second 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, through iterative orthogonalization, Hessenberg matrix update and least squares solution to the orthogonal basis matrix 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 Calculate updated approximate solution x i_vec , thereby obtaining a solution for the frequency-domain electromagnetic forward modeling calculation, and obtaining the electromagnetic response characteristics of the exploration object based on the solution for the frequency-domain electromagnetic forward modeling calculation. The present invention converts the left-hand matrix and right-hand terms of the large-scale sparse linear system obtained after discretization during the frequency-domain electromagnetic forward modeling calculation into a vectorized storage format. Through the concept of batch processing, sparse linear systems of multiple frequencies are solved simultaneously, and calculations are performed using vector instructions of modern CPUs. When using the generalized minimum residual method to solve the sparse linear system, each core operator (mainly including sparse matrix-vector multiplication, vector dot multiplication, and vector scalar multiplication) is expanded into a frequency-dependent vector calculation, and a convergence standard based on vectorized calculation is set to achieve the goal of solving sparse linear systems of multiple frequencies at one time. The present invention combines the computational characteristics and memory access characteristics of multi-frequency sparse linear systems, fully utilizes the long vector registers of the CPU, reduces the number of matrix memory accesses, solves the problems of reduced computational efficiency and low computational resource utilization in the existing technology, realizes fast and efficient multi-frequency frequency-domain electromagnetic forward modeling, and significantly improves the efficiency of frequency-domain electromagnetic forward modeling.

[0075] 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 calculations, such as parallel computing of matrix rows or calculations of other programs, which can save more computing resources.

[0076] During metal ore exploration, metal ore bodies often exhibit significant electrical anomalies (such as low or high resistance). However, their depth, shape, and size are complex and variable, making traditional forward modeling methods time-consuming and inadequate for large-scale exploration. The vectorized parallel solution method for frequency-domain electromagnetic forward modeling proposed in this paper significantly improves data processing speed, enabling rapid identification of electromagnetic response characteristics such as ore body location and electrical parameters. This method shortens the data processing cycle by over 50% while significantly reducing hardware costs. It also improves the accuracy of ore body boundary identification, reduces drilling blindness, and lowers exploration costs.

[0077] During oil and gas exploration, reservoirs often exhibit high-resistance anomalies. However, their electromagnetic responses are affected by complex stratigraphic structures (such as salt domes and faults), making it difficult for traditional forward modeling methods to rapidly simulate electromagnetic response characteristics such as multi-frequency and multi-scale electromagnetic field distributions. The vectorized parallel solution method for frequency-domain electromagnetic forward modeling proposed in this paper reduces reliance on high-performance computing clusters, efficiently simulates electromagnetic responses under complex geological conditions, shortens data processing cycles by over 50%, and significantly reduces hardware costs. It also improves the accuracy of reservoir boundary identification and reduces exploration risks. BRIEF DESCRIPTION OF THE DRAWINGS

[0078] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0079] Figure 1 1 is a schematic diagram of a vectorized parallel solution method for frequency domain electromagnetic forward modeling provided by an embodiment of the present invention;

[0080] Figure 2 1 is a schematic diagram of a vectorized parallel solution device for frequency domain electromagnetic forward modeling provided by an embodiment of the present invention;

[0081] Figure 3 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 DESCRIPTION

[0082] The present invention can be explained in detail through 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", "back", "left", "right", etc. to indicate directions or positional relationships, they only correspond to the drawings of this application for the convenience of describing the present invention, and do not indicate or imply that the device or element referred to must have a specific direction.

[0083] Example 1

[0084] refer to Figure 1 This embodiment provides a vectorized parallel solution method for frequency domain electromagnetic forward modeling, including the following steps:

[0085] S1, according to the vector register length and calculation accuracy of the CPU processor, the sparse linear system A of different frequencies is converted into i x i = b iConvert to vectorized calculation format, including sparse matrix A i The real matrix Re and sparse matrix A i The imaginary matrix Im i , reference complex matrix Im, frequency value f i , frequency vector f vec and the right end b i_vec Among them, A i is a sparse matrix of different frequencies, x i is the vector to be solved at different frequencies, b i is the right-hand term vector of different frequencies, i 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;

[0086] The sparse linear system A x = b Obtained by discretization of frequency domain electromagnetic forward calculation through finite element method or finite difference method;

[0087] A large-scale sparse linear system obtained by discretization of the finite element method (FEM) or finite difference method (FDM) in the same frequency domain electromagnetic forward calculation: 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, magnetic permeability, etc.) and discrete grid information; x is the vector to be solved, representing the value of the electromagnetic response characteristics to be solved (such as electric field or magnetic field) of the exploration object at the discrete grid nodes or units; b is the right-hand side vector, generated by source terms (such as current source or magnetic field source) and boundary conditions, representing external excitation or known conditions;

[0088] Different frequencies i The sparse linear system is represented by A i x i = b i ( i =1,2,……n) to represent the sparse linear system A x = b Transformed into a sparse linear system A with different frequencies i i x i = b i , where A i is a sparse matrix of different frequencies, A ifor nrow × nrow The complex matrix of ; x i is the vector to be solved at different frequencies, with dimensions of nrow * 1 vector; b i is the right-hand term vector of different frequencies, with dimension nrow * 1 Vector ；i is the serial number 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.

[0089] S11. Obtain the vector register length of the CPU processor;

[0090] The vector length of current mainstream CPU processors is usually 256 or 512 bits, and most of the newly released CPUs in recent years are 512 bits.

[0091] Generally, the vector length is indicated on the CPU processor model, that is, the vector length of the CPU processor, or by checking the CPU instruction set. For example, under the X86 architecture, if AVX256 is supported, the vector length is greater than or equal to 256; if AVX512 is supported, the vector length is greater than or equal to 512. This is existing technology and will not be repeated here.

[0092] S12. Calculate the number of vectors that can be simultaneously calculated by one CPU machine instruction based on the vector register length and calculation accuracy;

[0093] Assuming that the calculation precision used is 64-bit floating point and the vector register length of the computing platform is 512 bits, the number of vectors that can be calculated simultaneously by one machine instruction is 512÷64=8.

[0094] 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 is a sparse matrix A i Im represents the base complex matrix, f i is the frequency value;

[0095] Same grid with different frequencies i The generated sparse matrix A iIt has the following characteristics: A i Is a complex matrix, which can be expressed as A i =(Re+Im i ), where Re, Im i represent real matrices and imaginary matrices respectively.

[0096] Matrix A of different frequencies i The real part matrix Re is exactly the same; 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 Indicates frequency value, usually f i =0.1, 2.0, 1.0, 10.0, 100.0, etc.

[0097] S14, passing frequency value f i Constructing 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;

[0098] In step S12, it has been calculated that the number of vectorized instructions that can be calculated simultaneously for the current machine instruction is 8. Therefore, it is necessary to f vec The padded value is a multiple of the number of vectors. The padded value can be any value in the frequency value. The length of the padded vector is vec_num.

[0099] For example, the number of frequencies to be solved is 5. f i = 0.1, 2.0, 1.0, 10.0, 100.0, you need to f i Expand it into an array of length 8 and fill the tail, that is f vec =[0.1, 2.0, 1.0, 10.0, 100.0, f sup , f sup , f sup ]. f sup Indicates the fill value,f sup Usually the maximum value of the frequency values ​​involved in the calculation can be taken. f sup =100.0, vec_num = 8. If the number of frequencies to be solved is 10, then vec_num = 16, and so on.

[0100] 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 .

[0101] Through the right end vector b i Construct the right-hand side b i_vec , b i_vec for nrow A vector of rows, each element of which is a vec_num The long vector, nrow represents the total number of matrix rows, that is, the number of degrees of freedom to be solved.

[0102] use b i [ j ]( j =0,1,2,3……, nrow -1), indicating frequency i The corresponding right-hand term j Elements, reconstructed b i_vec become b 0 [ j ], b 1 [ j ], b 2 [ j ] continuous mode, similar to step S14, b i_vec Perform vector padding, and the padding value can be selected arbitrarily b i array.

[0103] For example, if the number of frequencies to be solved is 5, thenb i_vec The order of the elements of is as follows, each row represents a vector, where n = nrow -1, b sup [ j ] indicates a supplementary value.

[0104] {{ b 1 [0], b 2 [0], b 3 [0], b 4 [0], b 5 [0], b sup [0], b sup [0], b sup [0]};

[0105] { b 1 [1], b 2 [1], b 3 [1], b 4 [1], b 5 [1], b sup [1], b sup [1], b sup [1]};

[0106] …

[0107] { b 1 [n], b 2 [n], b 3 [n], b 4 [n], b 5 [n], b sup [n], b sup [n], b sup [n]}}.

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

[0109] S2. Initialization x i Approximate solution of x i_vec , and according to the maximum number of iterations max_iter and convergence tolerance tolerance Calculate residuals using vectorized calculation methods r i_vec and its second norm β vec Post-convergence update of approximate solution x i_vec ;

[0110] S21. Set the maximum number of iterations max_iter and convergence tolerance tolerance ;

[0111] Maximum number of iterations max_iter and convergence tolerance tolerance The value of is determined according to actual needs; preferably, the maximum number of iterations in this embodiment is max_iter Set to 100, the convergence tolerance tolerance Set to 10 -6 .

[0112] S22. Initialize the approximate solution x i_vec , x i_vec is a nrow A vector of rows, each element of which is a vec_num The long vector, x i_vec The value in can usually be 1;

[0113] For example, if the number of frequencies to be solved is 5, then x i_vec The order of elements is as follows, each row represents a vector, where n= nrow -1, x sup [ j ] indicates a supplementary value.

[0114] {{ x 1 [0],x 2 [0], x 3 [0], x 4 [0], x 5 [0], x sup [0], x sup [0], x sup [0]};

[0115] { x 1 [1], x 2 [1], x 3 [1], x 4 [1], x 5 [1], x sup [1], x sup [1], x sup [1]};

[0116] …

[0117] { x 1 [n], x 2 [n], x 3 [n], x 4 [n], x 5 [n], x sup [n], x sup [n], x sup [n]}}.

[0118] S23, vectorized calculation of initial residual r i_vec , r i_vec is a nrow A vector of rows, each element of which is a vec_num The long vector is calculated as follows;

[0119] r i_vec = bi_vec -(Re+ Im· f vec ) x i_vec

[0120] In the calculation process of sparse matrix multiplication by vector, the operation of multiplying each sparse matrix element by a vector is a multiplication of a vector and a vector;

[0121] For example, to calculate (Re + Im· f vec ) x i_vec The process:

[0122] (1) Take an i, j element of a matrix as an example;

[0123] (2) Expand Re[i][j] into a vec_num Vector Re vec [i][j];

[0124] (3) Calculate Im[i][j]· f vec Into a vec_num VectorIm vec [i][j];

[0125] (4) Get A vec [i][j]= Re vec [i][j]+ Im vec [i][j];

[0126] (5) Calculate A vec [i][j]· x i_vec [j,:]= v vec , v vec is a vec_num Vector.

[0127] The above calculations can all be performed by vectorized instructions of modern CPUs. From the calculation method, it can be seen that each matrix multiplication can complete all frequency-related matrix multiplication vector operations. The final result is r i_vec The order of elements is as follows, each row represents a vector, where n= nrow -1:

[0128] {{ r 1 [0], r 2 [0], r 3 [0],r 4 [0], r 5 [0], r 6 [0], r 7 [0], r 8 [0]};

[0129] { r 1 [1], r 2 [1], r 3 [1], r 4 [1], r 5 [1], r 6 [1], r 7 [1], r 8 [1]};

[0130] …

[0131] { r 1 [n], r 2 [n], r 3 [n], r 4 [n], r 5 [n], r 6 [n], r 7 [n], r 8 [n]}}.

[0132] S24. Calculation r i_vec The second norm of β vec , β vec is a dimension of vec_num vector;

[0133] β vec = sqrt ( r 0_vec 2 + r 1_vec2 + r 2_vec 2 +…… r n_vec 2 );

[0134] in r 0_vec represents a set of vectors,

[0135] { r 1 [0], r 2 [0], r 3 [0], r 4 [0], r 5 [0], r 6 [0], r 7 [0], r 8 [0]}.

[0136] The above calculations can be performed using vectorized computing instructions of modern CPUs.

[0137] The corresponding calculated β vec [i]= sqrt ( r i [0] 2 + r i [1] 2 + r i [2] 2 +…… r i [n] 2 ) (i=1, 2...8)

[0138] S25, convergence judgment, if β vec Every value in satisfies the convergence tolerance tolerance , then the output x i_vec .

[0139] 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 ;

[0140] S31. Initialize the orthogonal basis matrix V vec , the orthogonal basis matrix V vec Dimension nrow ×( max_iter +1) matrix, each matrix element is a vec_num A vector of length , used to store the orthogonal basis of the Krylov subspace, V vec The value of is initialized to 0;

[0141] S32, calculation V vec [:, 0], V vec [:, 0] = r i_vec / β vec ;

[0142] V vec The results stored in [:, 0] are as follows, each row represents a vector, where n= nrow -1:

[0143] {{ r 1 [0] / β [0], r 2 [0] / β [1], … r 8 [0] / β [7]};

[0144] { r 1 [1] / β [0], r 2 [1] / β [1], … r 8 [1] / β [7]};

[0145] …

[0146] { r 1 [n] / β [0], r 2 [n] / β [1], … r 8 [n] / β [7]}}.

[0147] In GMRES (Generalized Minimum Residual Method), the orthogonal basis matrix V vec It is a matrix composed of a set of orthogonal basis vectors generated by the Arnoldi process. V vec It is the auxiliary matrix constructed by GMRES to solve sparse linear systems.

[0148] S33. Initialize the Hessenberg matrix H vec , matrix H vec The construction dimension is max_iter ×( max_iter +1) matrix, each matrix element is a vec_num Vector of length, H vec The value of is initialized to 0.

[0149] The Hessenberg matrix H vec is the upper Hessenberg matrix;

[0150] S4, 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 Calculate updated approximate solution x i_vec , thus obtaining the solution of the frequency domain electromagnetic forward calculation.

[0151] Step S4 further includes: acquiring the electromagnetic response characteristics of the exploration object according to the solution of the frequency domain electromagnetic forward modeling calculation.

[0152] This step mainly performs the main loop iteration:

[0153] S41. Set the current number of iteration cycles iter =0;

[0154] S42. Calculate orthogonal vectors w i_vec= (Re + Im· f vec ) * V vec [: , iter ], w i_vec is a nrow A vector of rows, each element of which is a vec_num long vector;

[0155] The matrix multiplication vector operation involved can refer to the process of step S23 in step S2;

[0156] S43, correct orthogonalization, the specific steps are as follows:

[0157] S431, Settings k =0; Enter the loop;

[0158] S432, projection coefficient H vec [k][ iter ]calculate:

[0159] H vec [k][ iter ] = dot (( V vec [: , k ] ) T , w i_vec ), dot Represents the sum of two vector dot products;

[0160] S433, through w i_vec = w i_vec -H vec [ k ][ iter ] V vec [: , k ] to perform orthogonal correction;

[0161] S434, if k Less than iter+1 ,k = k + 1, Jump to S432, otherwise jump to S435;

[0162] S435, correction orthogonalization ends;

[0163] The matrix multiplication vector operation involved can refer to the process of step S23 in step S2;

[0164] S44. Calculate orthogonal vectors w i_vecThe second norm of and assign it to the Hessenberg matrix H vec ;

[0165] 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;

[0166] For detailed calculation method, refer to the vectorized two-norm calculation in step S24 in step S2.

[0167] S45. Normalized Hessenberg matrix H vec is an orthogonal basis matrix V vec ;

[0168] V vec [:, iter+1] = w i_vec / H vec [ iter +1][ iter ];

[0169] The process from S42 to S45 is the process of orthogonal basis expansion;

[0170] S46, through the least squares method to Hessenberg matrix H vec Solve to get the first vector array y i_vec and the second vector array e i_vec ;

[0171] This part cannot be vectorized and can only be solved sequentially; the solution process is as follows:

[0172] S461, Settings k =0; Enter the loop;

[0173] S462, Structure H k for( iter +2)×( iter +1) scalar matrix, H k Storage Hvec [: iter +2, : iter +1] in each vector element k +1 element;

[0174] For example, H0 stores H vec [: iter +2, : iter +1];

[0175] S463, structure e k for( iter +2) scalar array, initialize e k is 0; e k [0] = β [ k ];

[0176] S464, solve H by least squares method k y k =e k ,get y k is a ( iter +1);

[0177] S465, if k Less than vec_num,k = k +1 ， Jump to S462, otherwise jump to S466;

[0178] 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;

[0179] 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.

[0180] 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 convergence;

[0181] S471, calculation tmp vec =H vec [: iter +2, : iter +1] y i_vec -e i_vec ,

[0182] tmp vec is a iter +2 rows of vectors, each element of which is a vec_num Long vector.

[0183] S472, calculate Res vec , the residual value of the current iteration is as follows, Res vec is a vector of length vec_num:

[0184] Res vec = sqrt ( tmp 0_vec 2 + tmp 1_vec 2 + tmp 2_vec 2 +…… tmp n_vec 2 );

[0185] S473, judging the current residual Res vec Are all values ​​in less than the convergence tolerance tolerance If yes, go to step S49; if no, go to step S48;

[0186] S48, determine whether iter is less than the maximum number of iterations max_iter, if yes, jump to step S42, if not, jump to S49;

[0187] S49. Calculate the final solution and output the result.

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

[0189] Obtained x i_vec The solution in is the solution of the frequency domain electromagnetic forward modeling. x i_vec The storage information is as follows. x 1 [0]、 x 1 [1]…… x 1 [n] corresponds to a frequency of f 1 The solution is similar. x i [0]、 x i [1]…… x i [n] corresponds to a frequency of f i The value of the supplementary part can be discarded.

[0190] {{ x 1 [0], x 2 [0], x 3 [0], x 4 [0], x 5 [0], x sup [0], x sup [0], x sup [0]};

[0191] { x 1 [1], x 2 [1], x 3 [1], x 4 [1], x 5 [1], x sup [1], x sup [1], x sup[1]};

[0192] …

[0193] { x 1 [n], x 2 [n], x 3 [n], x 4 [n], x 5 [n], x sup [n], x sup [n], x sup [n]}}.

[0194] Finally, the electromagnetic response characteristics of the exploration object are obtained based on the solution of frequency domain electromagnetic forward modeling.

[0195] In existing technologies, magnetotellurics (MT) and controlled-source electromagnetics (CSEM) are both frequency-domain electromagnetic methods, widely used in metal exploration, oil and gas exploration, geothermal resource detection, and other fields. Frequency-domain electromagnetic forward modeling is one of the core technologies of frequency-domain electromagnetic methods. Traditional methods require solving for each frequency sequentially, which is inefficient; or they require multi-core parallel solutions, which requires additional computing resources.

[0196] During metal ore exploration, metal ore bodies often exhibit significant electrical anomalies (such as low or high resistivity). However, their depth, shape, and scale are complex and varied. Traditional forward modeling methods are time-consuming and difficult to meet the needs of large-scale exploration. The vectorized parallel solution method for frequency-domain electromagnetic forward modeling in this embodiment significantly improves data processing speed, enabling rapid identification of electromagnetic response characteristics such as ore body location and electrical parameters. This method shortens the data processing cycle by over 50% and significantly reduces hardware costs. It also improves the accuracy of ore body boundary identification, reduces drilling blindness, and reduces exploration costs.

[0197] During oil and gas exploration, reservoirs often exhibit high-resistance anomalies. However, their electromagnetic responses are affected by complex stratigraphic structures (such as salt domes and faults). Traditional forward modeling methods struggle to rapidly simulate electromagnetic response characteristics, such as multi-frequency and multi-scale electromagnetic field distributions. The vectorized parallel solution method for frequency-domain electromagnetic forward modeling in this embodiment reduces reliance on high-performance computing clusters, efficiently simulates electromagnetic responses under complex geological conditions, shortens data processing cycles by over 50%, and significantly reduces hardware costs. It also improves the accuracy of reservoir boundary identification and reduces exploration risks.

[0198] This embodiment can fully utilize the long vector register characteristics of modern CPU processors, can solve multiple frequency-related sparse linear systems at the same time, can greatly improve the calculation speed, and can significantly improve the efficiency of frequency domain electromagnetic forward modeling.

[0199] This embodiment proposes a vectorized parallel solution method for frequency domain electromagnetic forward modeling, including: S1, according to the vector register length and calculation accuracy of the CPU processor, the sparse linear system A of different frequencies is converted into i x i = b i Convert to vectorized calculation format; S2, initialization x i Approximate solution of x i_vec , and use vectorized calculation method to calculate the residual r i_vec and its second 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, through iterative orthogonalization, Hessenberg matrix update and least squares solution to the orthogonal basis matrix 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 Calculate updated approximate solution x i_vec. This embodiment converts the left-hand matrix and right-hand terms of the large-scale sparse linear system obtained after discretization in the frequency domain electromagnetic forward modeling calculation process into a vectorized storage format. Through the idea of ​​batch processing, the sparse linear systems of multiple frequencies are solved at the same time, and the calculations are performed through the vector instructions of modern CPUs. When using the generalized minimum residual method to solve the sparse linear system, each core operator (mainly including sparse matrix-vector multiplication, vector dot multiplication and vector scalar multiplication) is expanded into a frequency-related vector calculation, and a convergence standard based on vectorized calculation is set to achieve the goal of solving sparse linear systems of multiple frequencies at one time; this embodiment combines the computational characteristics and memory access characteristics of multi-frequency sparse linear systems, makes full use of the long vector registers of the CPU, reduces the number of memory accesses of the matrix, solves the problems of reduced computational efficiency and low computational resource utilization in the existing technology, realizes fast and efficient multi-frequency frequency domain electromagnetic forward modeling, and greatly improves the efficiency of frequency domain electromagnetic forward modeling.

[0200] 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 calculations, such as parallel computing of matrix rows or calculations of other programs, which can save more computing resources.

[0201] Example 2

[0202] refer to Figure 2 This embodiment provides a vectorized parallel solution device for frequency domain electromagnetic forward modeling, comprising the following units:

[0203] The vectorization unit is used to convert the sparse linear system A of different frequencies into i x i = b i Convert to vectorized calculation format, including sparse matrix A i The real matrix Re and sparse matrix A i The imaginary matrix Im i , reference complex matrix Im, frequency value f i , frequency vector f vec and the 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; bi is the right-hand term vector of different frequencies, with dimension nrow * 1 Vector ；i is the serial number 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;

[0204] The first initialization unit is used to initialize x i Approximate solution of x i_vec , and use vectorized calculation method to calculate the residual r i_vec and its second norm β vec Post-convergence update of approximate solution x i_vec ;

[0205] 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 ;

[0206] 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 Calculate updated approximate solution x i_vec , thus obtaining the solution of the frequency domain electromagnetic forward calculation.

[0207] Example 3

[0208] refer to Figure 3 , Figure 3FIG2 is a schematic diagram of the structure of a vectorized parallel solution device for frequency-domain electromagnetic forward modeling according to this embodiment. The vectorized parallel solution device 20 for frequency-domain electromagnetic forward modeling according to 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 of the aforementioned method embodiments are implemented. Alternatively, when the processor 21 executes the computer program, the functions of the modules / units in the aforementioned apparatus embodiments are implemented.

[0209] Exemplarily, the computer program can be divided into one or more modules / units, which are stored in the memory 22 and executed by the processor 21 to implement the present invention. The one or more modules / units can be a series of computer program instruction segments capable of performing specific functions, which are used to describe the execution process of the computer program in the vectorized parallel solver 20 for frequency-domain electromagnetic forward modeling. For example, the computer program can be divided into the modules described in Example 2. The specific functions of each module can be referred to the operating process of the device described in the above embodiment, and will not be repeated here.

[0210] The vectorized parallel solution device 20 for frequency-domain electromagnetic forward modeling may include, but is not limited to, a processor 21 and a memory 22. Those skilled in the art will appreciate that the schematic diagram is merely an example of the vectorized parallel solution device 20 for frequency-domain electromagnetic forward modeling and does not limit the vectorized parallel solution device 20 for frequency-domain electromagnetic forward modeling. The vectorized parallel solution device 20 may include more or fewer components than shown in the diagram, or may combine certain components or different components. For example, the vectorized parallel solution device 20 for frequency-domain electromagnetic forward modeling may further include input and output devices, network access devices, buses, and the like.

[0211] The processor 21 may be a central processing unit (CPU), other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA), other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or any conventional processor. The processor 21 serves as the control center of the vectorized parallel solver 20 for frequency-domain electromagnetic forward modeling, and utilizes various interfaces and circuits to connect various components of the vectorized parallel solver 20 for frequency-domain electromagnetic forward modeling.

[0212] The memory 22 can be used to store the computer programs and / or modules. The processor 21 implements the various functions of the vectorized parallel solver 20 for frequency-domain electromagnetic forward modeling by running or executing the computer programs and / or modules stored in the memory 22 and accessing the data stored in the memory 22. The memory 22 may primarily include a program storage area and a data storage area. The program storage area may store an operating system and at least one application required for a function (such as sound playback or image playback); the data storage area may store data generated based on the use of the mobile phone (such as audio data and a phone book). Furthermore, the memory 22 may include high-speed random access memory (RAM) and non-volatile memory, such as a hard disk, internal memory, a plug-in hard disk, a smart media card (SMC), a secure digital (SD) card, a flash card, at least one disk storage device, a flash memory device, or other volatile solid-state storage device.

[0213] If the modules / units integrated into the vectorized parallel solution device 20 for frequency-domain electromagnetic forward modeling are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the present invention can implement all or part of the processes in the above-mentioned method embodiments by instructing the relevant hardware through a computer program. The computer program can be stored in a computer-readable storage medium. When executed by the processor 21, the computer program can implement the steps of each of the above-mentioned method embodiments. The computer program includes computer program code, which can be in source code form, object code form, executable file, or some intermediate form. The computer-readable medium can include any entity or device capable of carrying the computer program code, a recording medium, a USB flash drive, a mobile hard drive, a magnetic disk, an optical disk, a computer memory, a read-only memory (ROM), a random access memory (RAM), an electrical carrier signal, a telecommunications signal, and a software distribution medium. It should be noted that the content contained in the computer-readable medium can be appropriately increased or decreased according to the requirements of legislation and patent practices in the jurisdiction. For example, in some jurisdictions, according to legislation and patent practices, computer-readable media do not include electrical carrier signals and telecommunication signals.

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

[0215] The parts of the present invention that are not described in detail are prior art. It is obvious to those skilled in the art that the present invention is not limited to the details of the above-mentioned exemplary embodiments, and that the present invention can be implemented in other specific forms without departing from the spirit or basic characteristics of the present invention. Therefore, no matter from which point of view, the embodiments should be regarded as exemplary and non-restrictive, and it is intended that all changes that fall within the meaning and scope of equivalent elements are included 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 of different frequencies is converted into i x i = b i Convert to vectorized calculation format, including sparse matrix A i The real matrix Re and sparse matrix A i The imaginary matrix Im i , reference complex matrix Im, frequency value f i , frequency vector f vec and the 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; b i is the right-hand term vector of different frequencies, with dimension nrow*1 Vector ；i is the serial number 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; 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 simultaneously calculated by one CPU machine instruction based on 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 is a sparse matrix A i Im represents the base complex matrix, f i is the frequency value; S14, passing frequency value f i Constructing 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 S2. Initialization x i Approximate solution of x i_vec , and use vectorized calculation method to calculate the residual r i_vec and its second 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 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 Calculate updated approximate solution x i_vec , thus obtaining the solution of the frequency domain electromagnetic forward calculation.

2. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 1 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 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 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 i_vec represents a set of vectors, β vec For a dimension vec_num vector, n= nrow -1.

3. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 2 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 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 , matrix H vec The construction dimension is max_iter ×( max_iter +1) matrix, each matrix element is a vec_num A vector of lengths.

4. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 3 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 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. Normalized Hessenberg matrix H vec is an orthogonal basis matrix V vec : V vec [:, iter+1] = w i_vec / H vec [ iter +1][ iter ]; S46, through the least squares method to 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; otherwise, jump to step S49. S49, according to the orthogonal basis matrix V vec and the first vector array y i_vec 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 .

5. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 4 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. Correction orthogonalization ends.

6. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 4, 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) scalar matrix, 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); 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, assembled 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, assembled 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.

7. The vectorized parallel solution method for frequency domain electromagnetic forward modeling according to claim 4, 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, judging the current residual Res vec Are all values ​​in less than the convergence tolerance tolerance If yes, go to step S49, if no, go to step S48.

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

9. 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 of different frequencies into i x i = b i Convert to vectorized calculation format, including sparse matrix A i The real matrix Re and sparse matrix A i The imaginary matrix Im i , reference complex matrix Im, frequency value f i , frequency vector f vec and the 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; b i is the right-hand term vector of different frequencies, with dimension nrow*1 Vector ；i is the serial number 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; specifically: S11. Obtain the vector register length of the CPU processor; S12. Calculate the number of vectors that can be simultaneously calculated by one CPU machine instruction based on 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 is a sparse matrix A i Im represents the base complex matrix, f i is the frequency value; S14, passing frequency value f i Constructing 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 The first initialization unit is used to initialize x i Approximate solution of x i_vec , and use vectorized calculation method to calculate the residual r i_vec and its second 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 Calculate updated approximate solution x i_vec .

Images