一种频率域电磁正演计算的矢量化并行求解方法及装置

1.一种频率域电磁正演计算的矢量化并行求解方法，其特征在于，具体包括以下步骤：S1、根据CPU处理器的矢量寄存器长度及计算精度将不同频率的稀疏线性系统A i x i = b i 转化为矢量化计算格式，包括稀疏矩阵A i 的实数矩阵Re、稀疏矩阵A i 的虚数矩阵Im i 、基准复数矩阵Im、频率值 f i 、频率矢量 f vec 及右端项 b i_vec ；其中，A i 为不同频率的稀疏矩阵，A i 为 nrow × nrow 的复数方阵； x i 为不同频率的待求解矢量，维度为 nrow*1 的矢量； b i 为不同频率的右端项矢量，维度为 nrow*1 的矢量 ；i 为不同频率的序号； i =1,2,……n，n为频率的个数，n为大于等于2的正整数； nrow 为需要求解的自由度的数量；步骤S1具体包括以下步骤：S11、获取CPU处理器的矢量寄存器长度；S12、根据矢量寄存器长度与计算精度计算CPU一条机器指令能够同时计算的矢量个数；S13、将不同频率的稀疏矩阵A i 转化为A i =(Re+Im i )，然后将Im i 转化为Im i = Im· f i ；其中，Re为稀疏矩阵A i 的实数矩阵，Im i 为稀疏矩阵A i 的虚数矩阵；Im表示基准复数矩阵， f i 为频率值；S14、通过频率值 f i 构建频率矢量 f vec ，并将 f vec 补齐至矢量个数的倍数，补齐后频率矢量 f vec 的矢量长度为vec_num；S15、通过右端项矢量 b i 构建右端项 b i_vec ， b i_vec 为 nrow 行的矢量，并将 b i_vec 每行补齐至矢量个数的倍数，补齐的值为右端项矢量 b i 中任一元素；S2、初始化 x i 的近似解 x i_vec ，并利用矢量化计算方法计算残差 r i_vec 及其二范数 β vec 后收敛更新近似解 x i_vec ；S3、初始化正交基矩阵 V vec 和Hessenberg矩阵H vec ，并计算 V vec [:, 0]， V vec [:, 0] = r i_vec / β vec ；S4、通过迭代正交化、Hessenberg矩阵更新及最小二乘求解对正交基矩阵 V vec 和Hessenberg矩阵H vec 进行更新并计算当前残差Res vec ，并在当前残差Res vec 收敛或迭代次数最大时根据正交基矩阵 V vec 计算更新近似解 x i_vec ，从而得到频率域电磁正演计算的解。 2.根据权利要求1所述的频率域电磁正演计算的矢量化并行求解方法，其特征在于，步骤S2包括以下步骤：S21、设置最大迭代次数 max_iter 及收敛容差 tolerance ；S22、初始化 x i 的近似解 x i_vec ， x i_vec 为一个 nrow 行的矢量，其中的每一个元素都是一个 vec_num 长的矢量；S23、矢量化计算初始残差 r i_vec ：r i_vec = b i_vec -(Re+ Im· f vec ) x i_vec ；r i_vec 为一个 nrow 行的矢量，其中的每一个元素都是一个 vec_num 长的矢量；S24、计算初始残差 r i_vec 的二范数 β vec ：β vec = sqrt ( r 0_vec 2 + r 1_vec 2 + r 2_vec 2 +…… r n_vec 2 )其中 r i_vec 表示一组矢量， β vec 为一个维度为 vec_num 的矢量，n= nrow -1。 3.根据权利要求2所述的频率域电磁正演计算的矢量化并行求解方法，其特征在于，步骤S3具体包括以下步骤：S31、初始化正交基矩阵 V vec ，所述正交基矩阵 V vec 为维度 nrow ×( max_iter +1)的矩阵，每一个矩阵元素上是一个 vec_num 长度的矢量；S32、计算 V vec [:, 0] ，V vec [:, 0] = r i_vec / β vec ；S33、初始化Hessenberg矩阵H vec ，矩阵H vec 构造维度为 max_iter ×( max_iter +1)的矩阵，每一个矩阵元素上是一个 vec_num 长度的矢量。 4.根据权利要求3所述的频率域电磁正演计算的矢量化并行求解方法，其特征在于，步骤S4具体包括以下步骤：S41、设置当前迭代循环次数 iter =0；S42、计算正交向量 w i_vec = (Re+ Im· f vec ) * V vec [: , iter ]， w i_vec 是一个 nrow 行的矢量，其中的每一个元素都是一个 vec_num 长的矢量；S43、修正正交化；S44、计算正交向量 w i_vec 的二范数，并赋值给Hessenberg矩阵H vec ：H vec [ iter +1][ iter ] = sqrt ( w 0_vec 2 + w 1_vec 2 + w 2_vec 2 +…… w n_vec 2 )，其中n= nrow -1；S45、标准化Hessenberg矩阵H vec 为正交基矩阵 V vec ：V vec [:, iter+1] = w i_vec / H vec [ iter +1][ iter ]；S46、通过最小二乘法对Hessenberg矩阵H vec 进行求解得到第一矢量数组 y i_vec 及第二矢量数组 e i_vec ；S47、根据Hessenberg矩阵H vec 、第一矢量数组 y i_vec 及第二矢量数组 e i_vec 计算当前残差Res vec 并检查是否收敛；S48、若当前残差Res vec 未收敛则判断iter是否小于最大迭代次数max_iter，若是则跳转步骤S42，若否则跳转步骤S49；S49、根据正交基矩阵 V vec 及第一矢量数组 y i_vec 对近似解 x i_vec 进行更新得到频率域电磁正演的解： x i_vec = x i_vec + V vec [:,iter+1] y i_vec 。 5.根据权利要求4所述的频率域电磁正演计算的矢量化并行求解方法，其特征在于，步骤S43具体包括以下步骤：S431、设置 k =0; 进入循环；S432、投影系数H vec [k][ iter ]计算：H vec [k][ iter ] =dot (( V vec [: , k ] ) T , w i_vec )， dot 表示两个矢量点乘求和；S433、通过 w i_vec = w i_vec -H vec [ k ][ iter ] V vec [: , k ]进行正交化修正；S434、如果 k 小于iter+1 ，k=k+1， 跳转至S432，否则跳转至S435；S435、修正正交化结束。 6.根据权利要求4所述的频率域电磁正演计算的矢量化并行求解方法，其特征在于，步骤S46具体包括以下步骤：S461、设置 k =0; 进入循环；S462、构造H k 为（ iter +2）×（ iter +1）的标量矩阵，H k 储存H vec [: iter +2, : iter +1]中每个矢量元素中的第 k +1个元素；S463、构造e k 为（ iter +2）的标量数组，初始化e k 为0；e k [0] = β [ k ];S464、通过最小二乘法求解H k y k =e k ，得到 y k 是一个（ iter +1）的标量数组；S465、如果 k 小于 vec_num，k=k +1 ， 跳转S462，否则跳转S466；S466、将 y k ，k=0，1，…… vec_num -1，组装成新的第一矢量数组 y i_vec ， y i_vec 是一个 iter +1行的矢量，其中的每一个元素都是一个 vec_num 长的矢量；S467、将 e k ，k=0，1，…… vec_num -1，组装成新的第二矢量数组 e i_vec ， e i_vec 是一个 iter +2行的矢量，其中的每一个元素都是一个 vec_num 长的矢量。 7.根据权利要求4所述的频率域电磁正演计算的矢量化并行求解方法，其特征在于，步骤S47具体包括以下步骤：S471、计算 tmp vec =H vec [: iter +2, : iter +1] y i_vec -e i_vec ，tmp vec 是一个 iter +2行的矢量，其中的每一个元素都是一个 vec_num 长的矢量；S472、计算Res vec ，当前迭代的残差值如下，Res vec 是一个vec_num长的矢量：Res vec = sqrt ( tmp 0_vec 2 + tmp 1_vec 2 + tmp 2_vec 2 +…… tmp n_vec 2 )；S473、判断当前残差Res vec 中的所有值是否均小于收敛容差 tolerance ，若是则跳转步骤S49，若否则跳转步骤S48。 8.根据权利要求1所述的频率域电磁正演计算的矢量化并行求解方法，其特征在于，稀疏线性系统A x = b 由频率域电磁正演计算通过有限元法或有限差分法离散后得到。 9.一种频率域电磁正演计算的矢量化并行求解装置，其特征在于，包括以下单元：矢量化单元，用于根据CPU处理器的矢量寄存器长度及计算精度将不同频率的稀疏线性系统A i x i = b i 转化为矢量化计算格式，包括稀疏矩阵A i 的实数矩阵Re、稀疏矩阵A i 的虚数矩阵Im i 、基准复数矩阵Im、频率值 f i 、频率矢量 f vec 及右端项 b i_vec ；其中，A i 为不同频率的稀疏矩阵，A i 为 nrow × nrow 的复数方阵； x i 为不同频率的待求解矢量，维度为 nrow*1 的矢量； b i 为不同频率的右端项矢量，维度为 nrow*1 的矢量 ；i 为不同频率的序号； i =1,2,……n，n为频率的个数，n为大于等于2的正整数； nrow 为需要求解的自由度的数量；具体包括：S11、获取CPU处理器的矢量寄存器长度；S12、根据矢量寄存器长度与计算精度计算CPU一条机器指令能够同时计算的矢量个数；S13、将不同频率的稀疏矩阵A i 转化为A i =(Re+Im i )，然后将Im i 转化为Im i = Im· f i ；其中，Re为稀疏矩阵A i 的实数矩阵，Im i 为稀疏矩阵A i 的虚数矩阵；Im表示基准复数矩阵， f i 为频率值；S14、通过频率值 f i 构建频率矢量 f vec ，并将 f vec 补齐至矢量个数的倍数，补齐后频率矢量 f vec 的矢量长度为vec_num；S15、通过右端项矢量 b i 构建右端项 b i_vec ， b i_vec 为 nrow 行的矢量，并将 b i_vec 每行补齐至矢量个数的倍数，补齐的值为右端项矢量 b i 中任一元素；第一初始化单元，用于初始化 x i 的近似解 x i_vec ，并利用矢量化计算方法计算残差 r i_vec 及其二范数 β vec 后收敛更新近似解 x i_vec ；第二初始化单元，用于初始化正交基矩阵 V vec 和Hessenberg矩阵H vec ，并计算 V vec [:,0]， V vec [:, 0] = r i_vec / β vec ；迭代计算单元，用于通过迭代正交化、Hessenberg矩阵更新及最小二乘求解对正交基矩阵 V vec 和Hessenberg矩阵H vec 进行更新并计算当前残差Res vec ，并在当前残差Res vec 收敛或迭代次数最大时根据正交基矩阵 V vec 计算更新近似解 x i_vec 。