一、概述

三对角线性方程组的求解是许多科学和工程计算中最重要也是最基本的问题之一。在核物理、流体力学、油藏工程、石油地震数据处理及数值天气预报等许多领域的大规模科学工程和数值处理中都会遇到三对角系统的求解问题。很多三对角线性方程组的算法可以直接推广到求解块三对角及带状线性方程组。由于在理论和实际应用上的重要性,近20年来三对角方程组的并行算法研究十分活跃。

大规模科学计算需要高性能的并行计算机。随着软硬件技术的发展,高性能的并行计算机日新月异。现今,SMP可构成每秒几十亿次运算的系统,PVP和COW可构成每秒几百亿次运算的系统,而MPP和DSM可构成每秒万亿次运算或更高的系统。

高性能并行计算机只是给大型科学计算提供了计算工具。如何发挥并行计算机的潜在性能和对三对角系统进行有效求解,其关键在于抓住并行计算的特点进行并行算法的研究和程序的设计与实现。另外,对处理机个数较多的并行计算系统,在设计并行算法时必须解决算法的可扩展性,并对可扩展性进行研究和分析。

二、问题的提出

设三对角线性方程组为

AX=Y (1)

式中:A∈Rn×n非奇异,αij=0, 。X=(x1,x2,…xn)T Y=(y1,y2,…yn)T。

此系统在许多算法中被提出,因此研究其高性能并行算法是很有理论和实际意义的。

三、并行求解三对角系统的直接解法

关于三对角线性方程组的直接求解已经有大量并行算法,其中Wang的分裂法是最早针对实际硬件环境,基于分治策略提出的并行算法。它不仅通信结构简单,容易推广到一般带状线性方程组的并行求解,而且为相继出现的许多其它并行算法提供了可行的局部分解策略。

近20年来求解三对角方程组的并行算法都是基于分治策略,即通过将三对角方程组分解成P个小规模问题,求解这P个小规模问题,再将这些解结合起来得到原三对角方程组的解。一般求解三对角方程组的分治方法的计算过程可分为3个阶段:一是消去,每台处理机对子系统消元;二是求解缩减系统(需要通信);三是回代,将缩减系统的解回代到每个子系统,求出最终结果。具体可分为以下几类:

(一)递推耦合算法(Recursive Doubling)

由Stone于1975年提出,算法巧妙地把LU分解方法的时序性很强的递推计算转化为递推倍增并行计算。D.J.Evans对此方法做了大量研究。P.Dubois和G.Rodrigue的研究表明Stone算法是不稳定的。

(二)循环约化方法(Cyclic Reduction)

循环约化方法由Hockey和G.Golub在1965年提出,其基本思想是每次迭代将偶数编号方程中的奇变量消去,只剩下偶变量,问题转变成求解仅由偶变量组成的规模减半的新三对角方程组。求解该新方程组,得到所有的偶变量后,再回代求解所有的奇变量。即约化和回代过程。由于其基本的算术操作可以向量化,适合于向量机。此方法有大量学者进行研究,提出了许多改进的方法。例如,Heller针对最后几步的短向量操作提出了不完全循环约化方法;R.Reulter结合IBM3090VF向量机的特点提出了局部循环约化法;P.Amodio针对分布式系统的特点改进了循环约化方法;最近针对此方法又提出对三对角方程组进行更大约化步的交替迭代策略。

(三)基于矩阵乘分解算法

将系数矩阵A分解成A=FT,方程Ax=b化为Fy=b和Tx=y两个方程组的并行求解。这种算法又可以分为两类:

1.重叠分解。如Wang的分裂法及其改进算法就属于这一类。P.Amodio在1993年对这类算法进行了很好的总结,用本地LU、本地LUD和本地循环约化法求解,并在1995年提出基于矩阵乘分解的并行QR算法。H.Michielse和A.Van der Vorst改变Wang算法的消元次序,提出了通信量减少的算法。李晓梅等将H.Michielse和A.Van der Vorst算法中的通信模式从单向串行改为双向并行,提出DPP算法,是目前最好的三对角方程组分布式算法之一。2000年骆志刚等中依据DPP算法,利用计算与通信重叠技术,减少处理机空闲时间取得了更好的并行效果。此类算法要求解P-1阶缩减系统。

2.不重叠分解。例如Lawrie Sameh算法、Johsoon算法、Baron算法、Chawla在1991年提出的WZ分解算法以及Mattor在1995年提出的算法都属于这一类。此类算法要求解2P-2阶缩减系统。

(四)基于矩阵和分解算法

将系数矩阵分解成A=Ao+△A,这类算法的共同特点是利用Sherman Morrison公式将和的逆化为子矩阵逆的和。按矩阵分解方法,这种算法又可分为两类:

1.重叠分解。这类算法首先由Mehrmann在1990年提出,通过选择好的分解在计算过程中保持原方程组系数矩阵的结构特性,具有好的数值稳定性,需要求解P-1阶缩减系统。

2.不重叠分解。Sun等在1992年提出的并行划分LU算法PPT算法和并行对角占优算法PDD算法均属于这一类。需要求解2P-2阶缩减系统。其中PDD算法的通讯时间不随处理机的变化而变化,具有很好的可扩展性。X.H.Sun和W.Zhang在2002年提出了两层混合并行方法PTH ,其基本思想是在PDD中嵌入一个内层三对角解法以形成一个两层的并行,基本算法是PDD,三对角系统首先基于PDD分解。PTH算法也具有很好的可扩展性。

四、并行求解三对角系统的迭代解法

当稀疏线性方程组的系数矩阵不规则时,直接法在求解过程中会带来大量非零元素,增加了计算量、通信量和存储量,并且直接法不易并行,不能满足求解大规模问题的需要。因此通常使用迭代法来求解一般系数线性方程组和含零元素较多三对角线性方程组。迭代法包括古典迭代法和Krylov子空间迭代法。

古典迭代法包括Jacobi、Gauss-Seidel、SOR、SSOR等方法。通常采用红黑排序、多色排序和多分裂等技术进行并行计算。

由于古典迭代法有收敛速度慢、并行效果不好等缺点,目前已较少用于直接求解大型稀疏线性方程组,而是作为预条件子和其它方法(如Krylov子空间方法)相结合使用。

Krylov子空间方法具有存储量小,计算量小且易于并行等优点,非常适合于并行求解大型稀疏线性方程组。结合预条件子的Krylov子空间迭代法是目前并行求解大型稀疏线性方程组的最主要方法。

给定初值X0,求解稀疏线性方程组AX=Y。设Km为维子空间,一般投影方法是从m维仿射子空间X0+Km中寻找近似解Xm使之满足Petrov-Galerkin条件

Y-AXm┻Lm

其中Lm为另一个维子空间。如果Km是Krylov子空间,则上述投影方法称为Krylov子空间方法。Krylov子空间Km(A,r0)定义为:

Km(A,r0)=span{r0,Ar0,A2r0,…,Am-1r0}

选取不同的Km和Lm就得到不同的Krylov子空间方法。主要算法包括四类:基于正交投影方法、基于正交化方法、基于双正交化方法、基于正规方程方法。

Krylov子空间迭代法的收敛速度依赖于系数矩阵特征值的分布,对于很多问题,直接使用迭代法的收敛速度特别慢,或者根本不收敛。因此使用预条件改变其收敛性,使中断问题可解,并加速收敛速度是需要的。目前人们研究的预条件技术可分为四类:采用基于矩阵分裂的古典迭代法作为预条件子、采用不完全LU分解作预条件子、基于系数矩阵近似逆的预条件子、结合实际问题用多重网格或区域分解作预条件子。对Krylov子空间和预条件Krylov子空间方法有详细的讨论。

预条件Krylov子空间方法的并行计算问题一直是研究热点,已提出了一系列好的并行算法。目前预条件Krylov子空间方法的计算量主要集中在矩阵向量乘上。虽然学者们做了大量的研究工作,但是还没找到效果好,又易于并行的预条件子。

需要特别指出的是,对于一般线性代数方程组的并行求解,其可扩展并行计算的研究已相对成熟,并已形成相应的并行软件库,如美国田纳西亚州立大学和橡树岭国家实验室研制的基于消息传递计算平台的可扩展线性代数程序库ScaLAPACK和得克萨斯大学开发的界面更加友好的并行线性代数库PLAPACK。我们借鉴其研究成果和研究方法,对三对角线性方程组并行算法的研究是有帮助的。

五、结语

三对角线性方程组的直接解法,算法丰富,程序较容易实现。但计算过程要增加计算量,并且大部分算法都对系数矩阵的要求比较高。迭代解法适合于非零元素较多的情况,特别是结合预条件子的Krylov子空间迭代法已成为当前研究的热点。

尽管三对角系统并行算法的研究取得了很多成果。但是还存在一些问题:直接法中,分治策略带来计算量和通信量的增加,如何减少计算量和通信量有待于进一步的研究;目前直接算法均基于分治策略,如何把其它并行算法设计技术,如平衡树和流水线等技术应用到三对角系统的并行求解中也是需要引起重视的方向;对于非对称系统还没找到一种通用的Krylov子空间方法;Krylov子空间方法的并行实现时仅考虑系数矩阵与向量乘,对其它问题考虑不够;以往设计的并行

算法缺乏对算法可扩展性的考虑和分析。

【参考文献】

[1]骆志刚,李晓梅,王正华.三对角线性方程组的一种有效分布式并行算法[J].计算机研究与发展,2000,(7)