并行算法设计的整体方案

上面已经对混凝土三维细观力学数值分析程序的结构进行了分析，对其中刚度矩阵装配算法与稀疏线性方程组的求解算法进行了改进，并对这些算法进行了详细的描述与分析。接下来将在该改进且优化过的程序的基础上，探讨对该问题的并行算法设计，包括整体设计方案，以及针对刚度矩阵装配、双门槛不完全Cholesky分解预条件、稀疏矩阵与向量相乘、稀疏向量相加等核心问题的高效并行算法。值得注意的是，这里所介绍的核心算法都是针对有限元问题的，所以对其他有限元问题，文中算法也具有借鉴意义。

在对与网格有关的问题进行并行计算时，区域分解是采用最多的方法之一。但对混凝土静动力学分析问题，在受外力位置以及加载时应力应变变化率较大的地方，在进行离散时，网格点应该比其他地方更加密集，所以在实际计算时，有限元网格点在物理空间上分布疏密不均，这使得在利用区域分解进行并行计算时，难以保证负载平衡，或为保证负载平衡，需要采用复杂的图分割技术（吴建平等，2004）。同时，如果需要对网格在有些位置进行加密，对负载平衡就会提出更严峻的挑战。为此，下面没有采用基于对网格进行简单剖分的方法来设计并行算法。

在进行并行算法设计时，整体上按有限元个数平均分配来进行任务划分，假设一共有p个处理器，将所有有限元分成p部分，每一部分中的有限元分配到一个处理器上。对一个给定的有限单元，如果分配到第k个处理器，则与之有关的计算在第k个处理器上进行。在对各个单元进行计算时，各个单元的计算之间是互不相关的，可以完全并行。

在各个单元对应的单元刚度矩阵形成后，需要利用这些信息形成整体刚度矩阵，在13.5.1节将详细描述该算法。在需要利用单元信息计算结点信息之处，先在各个处理器局部计算本地单元得出的结点信息，采用稀疏向量技术来存储，之后利用稀疏向量并行加法对各个处理器上的结点信息进行累加得到各个结点的整体信息，在13.5.4节将介绍稀疏向量相加的并行算法设计。

对维数很小的小数组，如材料参数矩阵AE（NPM，NM）、集中荷载矩阵FV（3，NL），及加载信息矩阵NF（NL）、损伤模型标识矩阵MNOD（NMID）、损伤模型参数矩阵PARAM（NM，5）、强化参数矩阵HCOEFF（NM，6）等，都需要从文件中输入，在并行算法设计时，从第0号进程输入，之后广播到所有进程，并在每个处理器上都保存一个副本。(www.xing528.com)

对与结点数有关的数组，主要也采用每个处理器上都保存一个副本的做法，如对输入的COOR坐标信息，在输入后直接广播到所有进程。在利用稀疏线性方程组求解结点位移时，与该线性方程组有关的右端向量与解向量是一种特殊情形。对这样的向量，直接以局部向量作为求解器的输入输出，当然矩阵也采用按行分块的方式进行分布。线性方程组求解时，采用并行预条件CG法（吴建平等，2004），其对应的基本串行算法参见13.3.1节和13.3.2节。由于对利用对角元比例化来增强对角优势的过程、向量更新、向量内积、与向量范数的并行算法设计十分直接，所以在这里不再详述，在13.5节中将详细介绍预条件构造过程、预条件作用过程、以及稀疏矩阵与向量相乘的并行算法设计。

在求出解向量之后，每个进程上只得到了局部分量，所以还需要进行收集操作，在每个进程上都得到整个解向量，以便利用结点位移计算各个单元的应变与应力增量。

采用如上所述的整体并行算法设计方案有两个优点：一是可以充分保证负载平衡，计算量比较大的部分要么对应于与有限单元有关的计算，要么对应于求解稀疏线性方程组，这两部分在设计方案中都几乎平均地在各个处理器上进行负载分配，所以负载平衡能得到保证；二是便于局部网格加密措施的采用。由于对网格加密后，可以再按单元分配与未知量分配相结合的方式设计并行算法，所以依然可以维持负载平衡。当然，由于未知量之间的数据相关性、以及有限单元之间的数据相关性都具有空间上的局部性，即只与网格中与之相邻的结点或有限单元有关，所以为减少通信时间，应该尽量保证空间上的相邻结点具有相邻的标号，同时使得空间上相邻的有限单元也具有相邻的标号。