一种用于计算矢量有限元方程组的不完全分解预处理方法

田 瑾，龚 利，史小卫，徐 乐

(1.西安电子科技大学天线与微波技术国家重点实验室，陕西西安710071;2.华东师范大学科技创新与发展战略研究中心，上海200062)

0 引言

矢量有限元(vector finite element method，FEM)是近40年来受到众多学者关注的一种求解边值问题的电磁计算方法，已被广泛应用于求解各类电磁问题［1－3］。金建铭［1］对该方法进行了较为全面、系统的介绍和说明，包括闭域和开域问题的几种求解方法。由于现实中研究的问题多是无限区域的，因此需要加入边界条件将开域问题转化为闭域问题来求解。通常采用的边界条件有3种:吸收边界条件(absorbing boundary condition，ABC)、理想匹配层(perfectly matched layer，PML)［2］和边界积分方程又叫边界元方法(boundary integral method，BI)［3］。前2种方法结合FEM方法离散目标得到的方程组的系数矩阵都是稀疏、对称的复方阵，而第3种方法结合FEM方法计算得到的方程组的系数矩阵是部分稀疏部分满秩的矩阵，针对后面的方法采用内观法分写方程组可以得到一个稠密系数复矩阵方程组和一个稀疏复矩阵方程组。

可以看出，不论采用何种边界条件，离散后的线性系统都可以归结为求解方程

(1)式中:A是稀疏/稠密、对称的复矩阵;x是待求列向量;b是列向量。通常采用直接法或迭代法求解(1)式。直接法得到的解在理论上是准确的，其缺点是高计算量和内存占用空间;迭代法可减少计算量提高计算速度，但若直接将迭代算法应用于大型有限元方程，收敛速度不能得到有效保证。由于迭代方法解的收敛速度主要由迭代矩阵谱特性决定，比较有效的求解方程组的方法是采用预处理方法改善迭代矩阵谱特性再结合迭代法进行计算［4－6］。文献［7］的数值结果表明，分解算法具有最快的收敛速度，也就是说，可以考虑采用不完全分解对方程组进行预处理，这样既能保证迭代法收敛速度，同时又不用付出完全分解所消耗的大量内存和时间。

直接采用分解法不能很好地利用缓存，而基于混合拓展乔里斯基(expanded Cholesky，ECM)［8］的多波前算法(multifrontal algorithm，MF)［9］(ECM/MF)从转化内存间接寻址为直接寻址的思想出发设计了一个向右算法，计算对称复系数矩阵方程组，该算法的分解操作是在2个稠密下三角矩阵内完成的，可以利用基本线性代数系统库函数(basic linear algebra subprograms，BLAS)［10］，达到很高的运算峰值。并且，在计算中，符号分解后引入缩放矩阵改善矩阵条件数以保证数值分解数值稳定性和之后的迭代计算的收敛速度。

预处理中采用的不完全分解法是在ECM/MF方法的基础上应用舍弃策略。文献［11］提出的舍弃策略在第j列运算中保留了分解因子的严格下三角部分的第j列中的coll个元素，coll为该列对应系数矩阵A严格下三角部分该列的非零元个数，这种算法可以预测存贮空间并考虑了元素数值的大小，但这种情况下可能舍弃元素过多，产生的预条件子显得过于粗糙，可能会导致迭代多次才能收敛或者不收敛。文献［12］提出的舍弃策略使分解因子拥有额外填充元，该算法保留了分解因子的严格下三角部分的第j列中的coll+p个元素，p是一个大于0的控制存贮空间的参数，然而这种情况下可能舍弃的元素又太少，从而使得预处理时间过长。本文结合上述2种方法的优点提出一种舍弃策略既能够保留部分填充元又能够控制预处理的时间。

针对(1)式中复系数矩阵对称的特点，本文考虑基于ECM/MF法构造一个高效率的不完全分解预处理方法。该算法是基于MF在超节点消去树的指导下在稠密矩阵内计算的原理，同时引入缩放矩阵以改善系数矩阵的条件数，并应用提出的舍弃策略生成预条件子，然后结合迭代法计算采用矢量有限元技术分析电磁问题产生的方程组。

1 预处理矩阵的构造

考虑电磁问题中(1)式的系数矩阵是复对称矩阵，A=LLT－E是一个不完全分解，L是下三角矩阵，则预处理线性系统可以表示为

预处理的目的是使得L－1AL－T收敛比原始矩阵A快，即迭代矩阵特征谱IN－L－1AL－T集中在0附近。下面首先简要说明ECM/MF分解算法，并给出舍弃策略，然后基于ECM/MF算法和舍弃策略提出不完全分解预处理方法并描述分解中用到的优化算法，最后给出用于迭代的方程组。

1.1 ECM/MF分解算法

ECM/MF分解算法在第k步计算的过程可由图1形象的给出，阴影部分为计算第k步时已被存贮在硬盘文件中的部分，黑色部分k为第k步正在分解的部分，黑色部分q为第k步正在更新的部分。采用ECM/MF分解法分解矩阵优势有2点:1)下三角矩阵L按列存贮，便于存贮和处理;2)每次只对一列元素进行分解操作，易于加入舍弃策略，从而构造不完全分解矩阵。

1.2 舍弃策略

为了保证预条件子不过于粗糙，同时避免预条件时间过长，本文以节点q为例提出的舍弃策略如下。

For j=1，…，m

dκ =τ‖low{A*j}‖F

判断Sq的第j列

判断1 如果该列有小于coll个非零元的

绝对值都≤dκ，则保留该列coll个最大非零元;

判断2 否则取t=min(coll2，coll+p)，

如果有大于t个非零元的绝对值都满足不小于某一

数值dκ，则保留该列t个最大非零元;

判断3 否则

保留该列所有≥dκ的非零元;

判断结束。

End for

其中j表示节点q相邻子节点的局部列编码;m是节点q相邻子节点数;* 表示A矩阵j列所有非零元的行号;Sq表示节点q的波前阵;coll2表示包括填充元在内的波前阵该列的非零元;p表示分解因子L中每列最大的填充元个数，一般取 p=nz(L)/N，nz(L)是分解因子L的非零元个数，N是矩阵K的阶数，K表示正在分解的矩阵;τ是舍弃容忍因子。预条件子的第j列保留的非零元素的个数至少为coll个元素，最多为t个元素。最后每列具体保留多少个，通过参数τ来确定，这样就不会浪费一些不必要的计算时间和存储空间。注意，如果选取的参数τ接近无限大的话，那么预条件子的第j列保留的元素个数就约等于coll个，这恰好是文献［11］提出来的算法。相反，如果选取的参数τ接近无限小或者为0，且coll2≥coll+p时，那么保留的元素个数就等于coll+p个，这恰好是文献［12］提出来的算法。

图1 ECM/MF分解法计算模式Fig.1 Decomposition model of ECM/MF

1.3 ECM/MF不完全分解算法

在ECM/MF算法基础上，引入缩放矩阵并应用提出的舍弃策略，构造出下列ECM/MF不完全分解算法。算法的流程如下。

1)根据系数矩阵的稀疏程度重排序方程组，采用反向卡特希尔－麦基(reverse Cuthill－Mckee，RCM)重排序稀疏矩阵系统，采用最小自由度(minimum degree，MMD)方法重排序中等稀疏和稠密矩阵系统，并令重排序后的矩阵为Aox=bo。

2)根据系数矩阵的正定程度计算对角矩阵P来缩放原系数矩阵，当原系数矩阵是正定系统时，P的对角元素为原系数矩阵的对角元素的幅值，否则，采用文献［13］的方法计算，使 ‖P－1AoP－1‖j∞∈［1 －μ，1+μ］，其中0μ1，‖‖j∞表示第j列的无穷范数，Ao表示A优化重排序后的矩阵。令缩放后的方程为 Ky=f，其中 K=P－1AoP－1，y=Px，f=P－1bo，bo表示b优化重排序后的向量。

3)采用ECM/MF方法计算方程组，以第j列为例，对波前阵F(j)用Submatrix－ECM分解，分解过程如下

(3)式中:Uj是分解进行中的ir×ir矩阵;Low{}表示取括号中矩阵包含对角元的下三角矩阵。图2给出了稠密矩阵空间Q1和Q2中的数据分解顺序。Q1和Q2内的箭头表示在左边节点波前阵计算结束后右边节点未更新的波前阵在该矩阵内替代左边节点波前阵数据，两个矩阵之间的箭头表示子节点更新父节点。箭头附近的数字表示计算操作顺序。节点1更新阵一旦形成立刻集成到其父节点4上，分解因子向量存入硬盘文件中，然后用节点2数据替换节点1数据，以此类推完成其他节点分解。也就是说，仅在2个稠密下三角矩阵内就可以完成所有节点分解操作。

图2 消去树和两个稠密矩阵中的分解顺序Fig.2 Decomposition order of elimination tree and Execution order of two dense matrices

4)对波前阵执行舍弃策略，则不完全分解后的L矩阵就是预条件子。

1.4 迭代求解

在预处理矩阵构造完毕后，采用迭代法求解矩阵方程Koz=fo，其中Ko=L－1KoL－1，z=Ly，fo=Lf。

基于上述缩放矩阵和舍弃策略的应用，本文实施的这种基于ECM/MF的不完全分解的预处理技术比传统的不完全分解预处理技术具有更好的稳定性和通用性。

2 数值算例

为了验证开域有限元方程组采用MF不完全分解预处理迭代法计算的正确性以及所编程序的可靠性，本文分别基于FEM/ABC，FEM/PML和FEM/BI方法采用提出的不完全分解预处理共轭梯度法(conjugate gradient，CG)计算了多种目标的双站或后向雷达散射截面图(radar cross section，RCS)。算法用Fortran编写，预处理CG法的收敛精度为1E－4，运行环境为 CPU 2.4 GHz，1.0 GByte RAM。

2.1 基于FEM/ABC方法计算RCS

分别计算尺寸为1.5λ×1.0λ×1.0λ理想导电矩形腔体后向散射和边长为0.755λ金属立方体的双站散射截面(θinc=180°，φinc=90°)问题。对于理想导电矩形腔体，虚拟边界为与腔体同中心，且半径r=1.35λ的球;对于金属立方体，虚拟边界为与腔体同中心，且半径r=0.855λ的球。图3给出理想导电腔体的后向散射截面(HH极化)和金属立方体双站雷达散射截面图，图3中实线和虚线是不完全分解预处理LG法的结果，空心方块是直接法求解的结果。由图3可以看出，本文方法与直接法计算结果一致。表1对比了本文不完全分解预处理CG法与直接法计算结果，表1中填充元个数比是不完全分解与完全分解的填充元个数比，可以看出，采用本文方法所需的内存远小于直接法，且计算时间明显优化。

图3 理想矩形导电腔体后向散射截面(HH极化)和金属立方体双站雷达散射截面图Fig.3 Backscatter RCS of a perfectly conducting rectangular inlet and bistatic RCS of a metallic cube

表1 不完全分解预处理CG法与直接法计算结果对比Tab.1 Comparison of incomplete factorization preconditioned CG and direct results

2.2 基于FEM/PML方法计算方向图

半波偶极天线PML模型如图4所示。图4中，PML层厚度为0.1λ，空气层厚度为0.15λ，方程组有237 824个未知量，2 283 402个非零元。图5给出了半波偶极天线的方向图，可见，采用不完全分解预处理CG法和直接求解的结果一致。本文方法与直接法的时间比值约为1/10，由此说明，在节省计算时间上，本文方法具有明显的优势。

2.3 基于FEM/BI方法计算RCS

半径r=0.5λ，介电常数εr=2.0的均匀介质球的双站雷达散射截面。采用内观法将FEM/BI方程组分写成两部分，稀疏矩阵部分采用RCM重排序，稠密矩阵部分采用 MMD重排序，矩阵乘用BLAS 2里面的DGMM命令计算。有限元部分有28 394个未知量，边界元部分有3 126个未知量。数值结果与Mie解析解比较见图6。

图6 均匀介质球的双站雷达散射截面图(θθ极化)Fig.6 Bistatic RCS of homogeneous dielectric sphere

3 结束语

本文提出了一种基于ECM/MF法的不完全分解预处理方法，该算法能够保证预条件子足够精细同时节约了构造时间。该算法改善了矩阵条件数和内存占用性能。数据结果表明，本文不完全分解预处理方法结合迭代法能够成功的应用于求解矢量有限元方程组，并在很大程度上改善存贮性能且提高计算速度。

［1］JIN J M.The finite element method in electromagnetics［M］.USA:Wiley－IEEE Press，1993.

［2］DAVIDSON D B，BOTHA M M.Evaluation of a spherical pml for vector FEM applications［J］.IEEE Transactions on Antennas and Propagation，2007，55(2):494－498.

［3］LI P.Coupling of finite element and boundary integral methods for electromagnetic scattering in a two－layered medium［J］.Journal of Computational Physics，2010，229(2):481－497.

［4］LI S S，RUI P L，CHEN R S.An effective sparse approximate inverse preconditioner for vector finite element analysis of 3d em problems［C］//Antennas and Propagation Society International Symposium［s.l.］:IEEE Press，2006.

［5］PING X W，CUI T.The Factorized sparse approximate inverse preconditioned conjugate gradient algorithm for finite element analysis of scattering problems［J］.Progress In Electromagnetics Research，2009(98):15－31.

［6］LI L，HUANG T，JING Y，et al.Application of the incomplete Cholesky factorization preconditioned Krylov subspace method to the vector finite element method for 3－D electromagnetic scattering problems［J］.Computer Physics Communications，2010，181(2):271－276.

［7］SHENG xin－qing，PENG zhen.Further cognition of hybrid FE/BI/MLFMA ［J］.Acta Electronica Sinica.2006，34(1):93－98.

［8］TIAN Jin，LV Zhi－qing，SHI Xiao－wei，et al.An efficient approach for multifrontal algorithm to solve non－positive－definite finite element equations in electromagnetic problems［J］.Progress In Electromagnetics Research，2009(95):121－133.

［9］QU Y，FISH J.Multifrontal incomplete factorization for indefinite and complex symmetric systems［J］.International journal for numerical methods in engineering，2002，53(6):1433－1459.

［10］LIU J W H.The multifrontal method for sparse matrix solution:Theory and practice［J］.Society for Industrial and Applied Mathematics，1992，34(1):82－109.

［11］JONES M T，PLASSMANN P E.An improved incomplete Cholesky factorization［J］.ACM Transactions on Mathematical Software(TOMS)，1995，21(1):5－17.

［12］LIN C J，MORÉ J J.Incomplete Cholesky factorizations with limited memory［J］.SIAM Journal on Scientific Computing，2000，21(1):24－45.

［13］QU Y，FISH J.Multifrontal incomplete factorization for indefinite and complex symmetric systems［J］.International journal for numerical methods in engineering，2002，53(6):1433－1459.