一种改进的基于线性有限元并行计算的追赶算法
开篇：润墨网以专业的文秘视角，为您筛选了一篇一种改进的基于线性有限元并行计算的追赶算法范文，如需获取更多写作素材，在线客服老师一对一协助。欢迎您的阅读与分享！
摘要 随着有限元计算规模的扩大，有限元并行计算的作用日益凸显。目前，对于有限元计算的并行处理的研究，主要集中在有限元计算的前置处理阶段，即有限元模型建立和网格划分阶段；研究以提高每个处理器所负责的子结构中的单元节点数的均匀程度及减少通信开销为目标。然而，鲜有对有限元分析的计算求解阶段的并行处理进行深入研究，而有限元并行计算绝大部分都在此阶段完成。本文就是针对该问题，对线性有限元并行计算求解阶段需要处理的刚度矩阵进行分析，依据刚度矩阵的特点和比较各种矩阵分解算法的优劣之后，提出一种优化的并行计算策略。经过验证，在对刚度矩阵的计算进行分布式并行处理后，有限元分析的计算速度有明显提高。
关键词 有限元并行计算；刚度矩阵；线性方程组；追赶算法；卷帘存储
中图分类号O241 文献标识码A 文章编号 1674-6708（2013）82-0112-04
1 概述
1.1 线性方程组系数矩阵分解及追赶算法
假设大型线性方程组用矩阵形式表示为，其中是系数矩阵。最为成熟的分解算法是高斯消元法。这种分解算法分为两个阶段。首先，用一些代数运算把化简成上三角形式的方程组，系数矩阵为，把这个方程组写成，其中是单位上三角矩阵。第二个阶段，使用回代法，对方程组进行求解。对这种方法的一种改进方法是，算法，就是将矩阵分解为上三角矩阵和下三角矩阵的分解算法，即。将存储在的下三角的位置，将存储在除主对角线以外的上三角的位置（因为的对角线元素都为1不用存储，默认为1）。高斯分解法适用于一般的稠密矩阵，对于正定矩阵，目前较为成熟的分解算法是Cholesky分解算法。它可以看成是分解法对对称正定矩阵的特例。即把一个矩阵分解成一个三角矩阵和它的转置的乘积，即，其中是上三角矩阵。显然，这种分解算法的计算量是Guess消元法的一半。对于三对角矩阵，Thomas在分解算法的基础上，提出了追赶算法。这种算法非常简单，计算量仅为次乘除运算，且追赶算法也是数值稳定的算法，是目前解三对角线性方程组的经典算法。
1.2 结构刚度矩阵
结构刚度矩阵的任一元素的物理意义是：结构中的第结点位移为单位值而其它结点位移皆为零时，需在第结点位移方向上施加的结点力的大小。与单元不同之处在于结构是单元的集合体，每个单元都对结构起一定的作用。由于单元刚度矩阵是对称和奇异的，由它们集成的结构刚度矩阵也是对称和奇异的，也就是说结构至少需给出能限制刚移的约束条件才能消除的奇异性，以便求得结点位移。
结构刚度矩阵是由单元刚度矩阵集合而成。因此，虽然总体单元数和结点数很多，结构刚度矩阵的阶数很高，但刚度系数中非零系数却很少，这就是刚度矩阵的大型和稀疏性。只要结点编号是合理的，这些稀疏的非零元素将集中在以主对角线为中心的一条带状区域内，即具有带状分布的特点。
综上所述，结构刚度矩阵的特点是具有奇异性、对称性、稀疏性和非零元素成带状分布的特点。如果单元节点的编号合理，可以确保结构刚度矩阵是正定三对角阵。
1.3 问题的提出
有限元分析是一种非常重要的数值分析方法，在工程技术及科学计算领域中得到了广泛的应用。然而，运用有限元分析方法对大型或超大型复杂结构进行有效的结构分析，所带来的计算量呈几何级数增长，通常方法是采用超级计算机进行计算。近年来，对有限元并行计算的研究受到研究者的关注，关注之一是改进有限元并行计算算法，使得这种大型复杂的结构分析在通常的分布式并行计算环境下，计算效率得到提高，适合于普通用户。
有限元分布式并行计算可分成三个阶段，前置处理、计算求解和后置处理。前置处理是建立有限元模型，完成单元网格划分；后置处理则是对处理结果进行分析，便于用户提取信息，了解计算结果。而大量的计算主要集中在求解部分，这一部分的计算本质上是大型线性方程组的求解，求解过程涉及线性方程组系数矩阵算法、追赶算法和结构刚度矩阵，求解的重点是如何对刚度矩阵进行分解。
本文将在现有线性方程组系数矩阵及追赶算法的基础上，结合结构刚度矩阵的特点提出一种有效刚度矩阵的分解策略，进而提出一种适合大型复杂结构分析的追赶算法，提高计算效率。
2.2 改进的适用于并行计算的追赶算法
2.2.1 现有的大型三对角矩阵的存储策略
1）块存储策略
假设是阶方阵，是结点机的台数。对于大型有限元分析，一般情况下，。假设。块存储（block storage）方案是，第一号结点机存储第一个行，第二号结点机处理第二个行，以此类推，最后一号结点机存储最后一个行，如果不是的整数倍，剩下的行都存在第一个节点中。由于矩阵靠前的行首先完成矩阵的分解，因此，每完成行，就增加一个空闲机器，这样出现了严重的负载不平衡现象。
2）单行卷帘存储策略
单行卷帘存储是将第行存放在号节点机中，比如第一号节点机会存储第l行，第+1行，第2+1行，以此类推。这样当第行分解完毕，第一个处理机才会停止运算。这样也最大程度的减少了负载不平衡的现象。但是，其通信开销也大大增加了。
3）多行卷帘存储策略
多行卷帘存储综合了两者的优势。这种算法是先将矩阵按行分块，每块包含多行。然后第块存储在第号结点机中。这样当第块分解完毕，第一个处理机才会停止运算。由于每块包含多行，通信开销要小于单行卷帘存储。
其中，它是除了最后一行外，其他元素均为的矩阵。由于是对角矩阵，可以分别把的最后一行除最后一个元素外均消为，记消去后的矩阵为，。从而，以及的最后一个元素构成一个新的小三对角阵方程组。对于这个小的三对角阵方程组可以在一个处理器上求解，然后将解发送给其他处理器，就可以求出愿问题的解。
2）算法描述
第一步：计算、，使；
第二步：对矩阵做变换；
第三步：计算形成小的三对角线性方程组，并求解之；
第四步：求解整个问题。
该算法具有良好的并行性，但是计算复杂度是串行算法的一倍。虽然增加了并行性，但是计算复杂度降低了算法的效率。因此，作者想在本算法的思路下进行改进，降低计算复杂度。
2.2.3 存储策略的改进
改进的多行卷帘存储（Improved Multiline Wrapped Interleaved Row Storage简称IMWIRS）
将数据存入第i个结点之后的存储情况。
第三步对于每个结点机，分别用追赶法计算小三对角阵和计算只含有一个元素的方阵，。其中，第一个节点不用处理，最后一个节点不用处理。
2）算法描述
在2.2.1和2.2.2描述的分解策略和映射技术的基础上，给出本文提出的算法。在存储上，在多行卷帘存储的基础上，节省内存开销，只存储对角线附近的元素。在分解策略和映射技术上，将大型三对角矩阵分解成小的方阵，每个处理器需要处理一个小的三对角矩阵和两个只含有一个元素的方阵。然后每个处理器将计算结果发送给主处理器，计算原始问题的结果。具体算法如下所述。
假设并行计算机的网络拓扑为主从架构，要计算的线性方程组用矩阵表示为。
第一步：在主计算机上输入要计算的大型刚度矩阵的阶数，假设我们的程序要将矩阵按行划分为块（或者说，有个从节点），那么主计算机按照多行卷帘存储策略计算出每块的行数。
第二步：输入要计算的刚度矩阵的元素和向量。在输入的时候，主计算将系数矩阵Ａ的第l至行的元素，向量分配给第一个从结点机，第到行的元素，向量分配给第二个从结点机，以此类推，最后剩下的不足行的元素和向量分配给第一个结点机。注意，在分配的时候，主计算机只处理三对角元素，即下标的元素。
第三步：在每个从计算机上开辟一个的数组，将主计算机分配的元素按照IMWIRS（见2.2.1节）方法存储。
第四步：对于每一个从结点机在逻辑上将数组中的元素分为三个阶方阵，数组的第一个元素看成的右下角元素，数组的最后一个元素看成是的左上角元素，数组的其他元素是中的元素。映射方法见
4 结论及展望
实验结果证明，本算法在提高计算效率、降低计算复杂度和节约内存上都有不错的效果。
1）本算法的内存占用为一般算法的3/n；
2）交互开销极小，计算复杂度远远小于经典的“分裂法”（分裂法介绍见2.2.2）；
3）随着处理器的增加，加速比增加，但是到了16台之后，加速比增加不明显。并行效率下降，是由于增加了处理器的台数，也增加了通信开销；
4）随着矩阵规模的扩大，并行效率增加。
本算法没有考虑局域网的体系结构的变化，以后可能会在这个方面展开研究。
参考文献
[1]AnanthGrama，AnshulGupta&GeorgeKarypis，Vipin Kumar等.Introduction to Parallel Computing[M].北京：机械工业出版社，2003.
[2]C.Xavier，S.S.IyengarIntroduction to Parallel Algorithms[M].北京：机械工业出版社，2004.
[3]Sun Wenyu，Du Qikui，Chen Jinru.计算方法[M].北京：科学出版社，2007.
[4]Wang Xingbo，ZhongZhihua.Generalized Thomas algorithm for solving cyclic tridiagonal equations[J].计算力学学报，2004（2）：73-76.
[5]Hee-Seok Kim，Shengzhao Wu，Li-wen Chang.A Scalable Tridiagonal Solver for GPUs[C].Parallel Processing （ICPP），2011 International Conference on.2011（9），444-453.
[6]李文强，刘晓.追赶法并行求解循环三对角线性方程组[N].科技导报，2009（5）：90-93.