典型結(jié)構(gòu)大型線性代數(shù)方程組的求解

1、文獻(xiàn)綜述1前言典型結(jié)構(gòu)大型線性代數(shù)方程組的求解是許多應(yīng)用領(lǐng)域的基礎(chǔ)，如：結(jié)構(gòu)分析、電子工程、油藏模擬、計(jì)算機(jī)輔助幾何設(shè)計(jì)、大氣污染研究、化學(xué)工程和經(jīng)濟(jì)模 型模擬、核物理和計(jì)算流體力學(xué)、數(shù)值天氣預(yù)報(bào)等。在數(shù)學(xué)物理及工程技術(shù)領(lǐng)域， 如微分方程的求解、多項(xiàng)式插值、網(wǎng)絡(luò)、系統(tǒng)控制等方面也常會(huì)碰到的大型、分 塊三對(duì)角矩陣為系數(shù)陣的線性方程組的求解問題。一般，動(dòng)態(tài)過程的數(shù)學(xué)模型由 偏微分方程描述，而偏微分方程的離散化通常導(dǎo)出大型線性方程組，它們可能是對(duì)稱或非對(duì)稱的大型稀疏線性方程組，也可能是結(jié)構(gòu)化的大型稀疏線性方程組。 甚至于對(duì)于依賴于時(shí)間的非線性問題，其全局計(jì)算的中間步驟也需要對(duì)線性方程 組的求解。長(zhǎng)期

2、以來，伴隨著計(jì)算環(huán)境的不斷變化，人們對(duì)于求解各類大型線性 方程組的適應(yīng)新的計(jì)算環(huán)境的新方法的探求從來也沒有停止過。目前，分布式存儲(chǔ)并行處理機(jī)系統(tǒng)己經(jīng)成為許多科學(xué)和工程問題的計(jì)算環(huán)境，成為求解重大挑戰(zhàn)性問題的首選工具；工作站機(jī)群（NOWs）和PC機(jī)群作為具有良好性價(jià)比的分布 式存儲(chǔ)并行處理機(jī)系統(tǒng)已廣泛應(yīng)用于各類科學(xué)和工程計(jì)算問題。典型結(jié)構(gòu)大型線性方程組的解法從總體上說可分為直接法和迭代法兩大類。 求解具有結(jié)構(gòu)化系數(shù)矩陣的大型線性方程組的研究近年主要集中在直接法，而迭代法近年來已發(fā)展為求解一般大型稀疏線性方程組的主要方法。本文所研究的內(nèi)容如下：考慮大型線性方程組Ax=b, AWRn x、bw Rn

3、 ,其中A為三對(duì)角或塊三對(duì)角系數(shù)矩陣，探討分布式存儲(chǔ)環(huán)境下求解大型線性方程組的高效并行算法。在科學(xué)與工程問題中經(jīng)常遇到的許多微分方程， 經(jīng)過適當(dāng)差分或有限元離散 而形成系數(shù)矩陣是塊三對(duì)角的線性方程組， 它們的求解是高性能并行計(jì)算的重要 課題之一。目前針對(duì)求解塊三對(duì)角線性方程組的并行算法的研究已經(jīng)有了一些成 果，通過對(duì)系數(shù)矩陣進(jìn)行分解與近似處理， 構(gòu)造了具有良好的并行性的算法。 借 助現(xiàn)有的并行工具環(huán)境，進(jìn)一步構(gòu)造出了并行效率更高的并行求解算法。2研究現(xiàn)狀求解典型結(jié)構(gòu)三對(duì)角線性代數(shù)方程組有多種方法，其解法總體可分為直接法 和迭代法兩大類。迭代法（ iterative methods） 1,2主要

4、包括 Jacobi迭代、 Gauss-Seidel迭代、逐次松弛迭代法（SOR）,直接法包括高斯消元和幾類并行算 法。迭代法Jacobi迭代因各個(gè)分量的修正相互獨(dú)立而具有十分明顯的內(nèi)在并行計(jì)算特 性。其主要優(yōu)點(diǎn)是方法簡(jiǎn)單，然而它并不常是收斂的，收斂時(shí)速度常較慢。在研 究如何提高收斂速度的基礎(chǔ)上，1983年，Missirlsi提出了并行Jacobi型方法，并 討論了它的收斂性。胡家贛等把它推廣到兩參數(shù)的情形，稱之為兩參數(shù)并行 Jacobi 型方法3。對(duì)于Gauss-Seidel迭代，因充分利用上次求出的新值，可加快收斂速度， 正因?yàn)槊看吻笾刀家玫缴洗蔚男轮?，使它不容易并行。?duì)于SOR迭代法來說

5、，由于各分量的計(jì)算是逐個(gè)相關(guān)的，因此，一般認(rèn)為 SOR迭代法不適合并行處理，其內(nèi)在并行性遠(yuǎn)不如Jacobi迭代。由于SOR多用于有限差分或有限元方法導(dǎo)致的大型稀疏方程組的求解。因此，利用系數(shù)矩陣零元素或非零元素的特殊分布，采用紅黑或多色排序成為實(shí)現(xiàn)SOR并行處理的有效途徑。然而，如何找到合適的“彩色模板”并保持自然排序下的收斂速度卻是 一個(gè)問題。蔡放等4對(duì)SOR方法通過改造提出的向量化 SOR算法，在一定條件 下具有較好的并行化計(jì)算性能。后來，呂全義在文獻(xiàn)5中對(duì)BSOR方法通過改進(jìn)引入加速因子和松弛因子，使收斂速度相同，但降低了迭代次數(shù)。接著，崔喜寧6在此基礎(chǔ)上，結(jié)合P&方法，將內(nèi)迭代

6、采用PH方法，使矩陣的分裂在一定程度上有大的改進(jìn)，使算法更具靈活性，并行效率也很高。以上幾種基本迭代方法是進(jìn)行并行迭代的基礎(chǔ)，充分了解其并行再借鑒用行算法進(jìn)行并行程序設(shè)計(jì)，在這些基礎(chǔ)上研究新的算法并重新獲得快速的收斂速 度。肖曼玉和呂全義在文獻(xiàn)7中，提出了一種基于 Galerkin原理求解塊三對(duì)角 線性方程組的Arnoldi并行算法，通過選取適當(dāng)?shù)淖涌臻g，使算法只在相鄰處理 機(jī)間有通信，因而具有很好的并行性，而且證明了該算法的收斂性。在HPrx2600 集群上進(jìn)行數(shù)值計(jì)算，結(jié)果表明，加速比呈線性增加，并行效率達(dá)到90%以上。直接法托馬斯(Thomas)算法8是一種對(duì)于三對(duì)角線性方程組的特殊的高

7、斯消元 方法，也是求解三對(duì)角線性方程組首先想到的解法。Thomas算法的求解過程可以概括為兩個(gè)步驟(removing和backward),首先從前往后依次消去對(duì)角線下方 的非零元素；再反向回代解出的值，從后往前依次解出所有的未知變量。 此算法 也容易擴(kuò)展到塊三對(duì)角矩陣9的應(yīng)用中。雖然托馬斯算法是申行計(jì)算機(jī)上的最 快算法，但由于后面的每一步都要依賴前一步的計(jì)算結(jié)果，它是不可并行的。三對(duì)角線性方程組和塊三對(duì)角線性方程組的并行算法，其研究始終非?；?躍?？偨Y(jié)以往算法，可將其歸為以下幾類：(1)循環(huán)遞減算法(cyclic reduction), (2)遞推倍增算法(recursive doubling

8、), (3)矩陣分解算法(partition method)。 雖然有Amodio10將CR改進(jìn)后用到超立方體結(jié)構(gòu)機(jī)器，總體而言，遞推倍增算 法和循環(huán)遞減算法(CR)是適應(yīng)向量計(jì)算機(jī)或共享內(nèi)存并行機(jī)的并行算法。Lambiotte and Voigt 1提出了基于CDC STAR-100計(jì)算機(jī)的循環(huán)遞減算法，按照 這種算法，經(jīng)過一次遞減后，原來的線性方程組中的奇下標(biāo)未知數(shù)都被消去了， 而留下了所有偶下標(biāo)的未知數(shù)變量，于是得到原有線性方程組的一半規(guī)模的同結(jié) 構(gòu)三對(duì)角線性方程組11。對(duì)新得到的減半線性方程組重復(fù)采用循環(huán)遞減算法， 直到最后剩下只含兩個(gè)未知數(shù)的線性方程組，并解出這兩個(gè)未知數(shù)。后續(xù)步驟就

9、 是，回代解出的變量至其上一步遞減的線性方程組，遞歸求解得到每個(gè)未知量的解值。如此，回代求解出整個(gè)三對(duì)角線性方程組的解值12 0近年來，隨著計(jì)算環(huán)境的發(fā)展，循環(huán)遞減算法被許多人應(yīng)用到GPU上實(shí)現(xiàn)13,14。Stone15,16方次提出了遞推倍增算法。在遞推倍增算法中，涉及到系數(shù)矩陣的 LU分解，以及 順向和逆向遞推方法。Wang17提出的分裂法屬于矩陣分解算法， 適用于分布存 儲(chǔ)環(huán)境，受到關(guān)注，Michielse和van der Vorst18改進(jìn)了 Wang的分裂法。遲利 華和李曉梅19在Michielse和van der Vorst算法18的基礎(chǔ)上，提出了雙向并行 分裂法(DPP算法)。另

10、外，Bondeli18提出了一種基于分治思想的算法；Mu11er 及Scheerer21fc 1991年提出了一種將三對(duì)角方程組用行算法并行化的一般方 法。至于塊三對(duì)角系統(tǒng)和帶狀系統(tǒng)方面，Rodrigue等22曾將奇偶約化法推廣到 帶狀系統(tǒng)的并行求解；Meier23將Wang的劃分算法17拓展到帶狀方程組； Kapur和Brown24提出了一種適用于可重構(gòu)陣列計(jì)算機(jī)的塊三對(duì)角線性方程組 并行算法；van der Vorst18于1987年基于不完全分解提出了一種向量并行機(jī)上 求解大型三對(duì)角和塊三對(duì)角方程組的方法；Ruggiers及Galligani25提出了一種基于迭代法和預(yù)條件子的塊三對(duì)角方

11、程組的并行解法。對(duì)于循環(huán)塊三對(duì)角線性方程組，文獻(xiàn)26中討論了適用于共享主存向量機(jī)的 并行算法；胡慶豐、何新芳、李曉梅1791提出了以分塊壓縮存儲(chǔ)形式直接求解的分塊追趕法及其在向量機(jī)上的并行計(jì)算方案。Chung等27給出了一種基于“分 治”思想的并行算法，可用于分布存儲(chǔ)計(jì)算環(huán)境。Toeplitz系統(tǒng)在數(shù)學(xué)、數(shù)字信號(hào)處理和ARMA模型中均有廣泛應(yīng)用，Toeplitz 三對(duì)角和Toeplitz循環(huán)三對(duì)角線性方程組的求解，是具有結(jié)構(gòu)化系數(shù)矩陣的大型 稀疏線性方程組求解中的重要課題。對(duì)于這類方程組，Evans提出了一種基于系數(shù)矩陣分解的快速算法28,是目前求解這類方程組的最快的用行算法；Buckley

12、給出了一種算法29,它是由通常的高斯消去法改進(jìn)而來的；關(guān)于 Toeplitz循環(huán) 三對(duì)角線性方程組，趙自春、李曉梅提出了一種適用于向量并行機(jī)的并行算法 30 o關(guān)于Toeplitz三對(duì)角線性方程組，趙自春、李曉梅提出了兩種向量并行計(jì) 算機(jī)上適用的并行算法31 o Evans和Yousif在文獻(xiàn)32、文獻(xiàn)33基礎(chǔ)上提出了 適用于共享存儲(chǔ)多處理機(jī)的一種并行算法，它是循環(huán)奇偶約化法變化和改進(jìn)而得 的。成禮智和蔣增榮34對(duì)帶狀（塊）Toeplitz方程組進(jìn)行討論，給出了向量機(jī) 上的一個(gè)快速并行算法。Poisson方程的求解，導(dǎo)致一類特殊的實(shí)對(duì)稱塊三對(duì)角線性方程組的求解， 這一問題也曾受到廣泛關(guān)注，它的

13、串行算法和向量或共享共享存儲(chǔ)計(jì)算機(jī)上的并 行算法得到了充分的研究。Buzbee等35討論了求解這類方程組的直接解法：矩陣分解法（MD）和塊循環(huán)遞減法（BCR）,以及結(jié)合MD與BCR的FACR （l）, Buzbee等35在同一文獻(xiàn)中還對(duì)MD方法和BCR方法應(yīng)用于Poisson方程的情 形進(jìn)行了全面的討論，涉及到基于FFT的MD方法以及具有數(shù)值穩(wěn)定性的Buneman算法36。Sweet37推廣Buneman的算法，提出了 n為一般正整數(shù)時(shí) 的循環(huán)遞減法。Swarztrauber38提出了求解Poisson方程的近似循環(huán)遞減（ACR） 方法。在共享存儲(chǔ)環(huán)境，上述方法都易于并行實(shí)現(xiàn)，文獻(xiàn)26對(duì)基于F

14、FT的MD方法和Buneman算法的并行實(shí)現(xiàn)進(jìn)行了討論。三角形矩陣是一種特殊結(jié)構(gòu)的矩陣。高效率地并行求解三角形方程組具有非 常重要的意義，這是因?yàn)橄禂?shù)矩陣 A為一般稠密矩陣的大型線性方程組的各種 直接解法大都基于化系數(shù)矩陣為三角形矩陣的思路來處理。分布式環(huán)境下求解三 角形方程已有研究者進(jìn)行過許多有益的探討，如 Heath和Romine39, Eisenstat 等40, Li 和 Coleman41,42, Fiebach43等。Li 和 Coleman于 1988 年提出了 一種求解系數(shù)矩陣以列（或行）卷簾方式分布時(shí)的三角形方程組并行解法41,1989年他們又對(duì)算法進(jìn)行了改進(jìn)42 0 Fie

15、bach于1996年提出了一種適于網(wǎng)絡(luò)拓 撲為二維網(wǎng)格的分布存儲(chǔ)多計(jì)算機(jī)系統(tǒng)上求解三角形方程組的循環(huán)塊算法43 o3總結(jié)許多計(jì)算問題都需要求解三對(duì)角線性方程組，這方面最為普通的例子當(dāng)屬橢圓型偏微分方程的差分格式求解了。三對(duì)角線性方程組的并行算法是數(shù)值并行算 法研究中最重要的問題之一，其研究始終非?；钴S。伴隨著計(jì)算機(jī)體系結(jié)構(gòu)的發(fā) 展，人們不斷地提出和改進(jìn)了許多三對(duì)角線性方程組的并行算法，直接法中由 Wang17提出的分裂法是貫徹分而治之原則的成功例子，受到廣泛重視。Michielse和van der Vorst18改進(jìn)了 Wang的分裂法，減少了通信開銷，提高了 計(jì)算與通信的重疊程度。遞M倍增算

16、法( recursive doubling)和循環(huán)遞減算法(cyclic reduction)同樣是備受關(guān)注的典型并行方法，已經(jīng)被多種并行計(jì)算技術(shù) 實(shí)現(xiàn)。塊(循環(huán))三對(duì)角線性方程組的求解同樣也是諸多應(yīng)用問題的重要組成部 分。例如，周期的樣條插值就導(dǎo)致對(duì)角占優(yōu)的循環(huán)三對(duì)角線性方程組的求解，當(dāng)邊界條件為周期邊界條件時(shí)，一些偏微分方程的離散化也可能導(dǎo)致循環(huán)三對(duì)角線 性方程組的求解。4參考文獻(xiàn)1. J. J. Lambiotte Jr., R. G. Voigt. The solution of tridiagonal linear systems on the CDC STAR-100 comput

17、er. ACM Trans. Math. Softw. 1,308329 (1975)2. Huy Nguyen, Matthew A. Beauregard, Ron Morgan. Improving the Speed of Convergence of GMRES for Certain Perturbed Tridiagonal Systems. In: IEEE 45th Southeastern Symposium on System Theory(SSST), pp. 667 (2013)3. 張寶琳，谷同祥，莫?jiǎng)t堯.數(shù)值并行計(jì)算原理與方法.北京:國(guó)防工業(yè)出版 社，19994.

18、 蔡放，方透，李彬.具有完全并行度的SOR算法.高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào)， 1:11 T4 (2002)5. 呂全義，葉天麒.系數(shù)矩陣為塊三對(duì)角的線性方程組的并行算法.西北工業(yè)大學(xué)學(xué)報(bào),14(2): 314 W18 (1996)6. Xinning Cui, Qunayi Lu. A Parallel Algorithm for Block-tridigaonal Linear Systems, Applied Mathematics and Computation, 2: 11071114 (2006)7. 肖曼玉，呂全義.塊三對(duì)角線性方程組的一種有效并行算法J.計(jì)算機(jī)應(yīng) 用與軟件,23(6):

19、 107-108 (2006)8. L.H. Thomas. Elliptic problems in linear difference equations over a network.Watson Scientific Computing Laboratory Report, Columbia University (1949)9. S.P. Hirshman, K.S. Perumalla, V.E. Lynch, R. Sanchez. BCYCLIC: a parallel block tridiagonal matrix cyclic solver. J. Comput. Phy

20、s. 229, 6392-6404 (2010)10. AmodioP, MastronardiN. A parallel version of the cyclic reduction algorithm on a hypercube, Parallel Computing, 19: 127 31281 (1993)11. Tomohiro Sogabe, Moawwad El-Mikkawy. Fast block diagonalization of k-tridiagonal matrices. Applied Mathematics and Computation 218, 2740

21、 -2743 (2011)12. Sudip K. Seal, Kalyan S. Perumalla, Steven P. Hirshman. Revisiting parallel cyclic reduction and parallel prefix-based algorithms for block tridiagonal systems of equations. J. Parallel Distrib. Comput. 73, 273280 (2013)13. Yao Zhang, Jonathan Cohen, John D. Owens. Fast tridiagonal

22、solvers on theGPU. In: Proceedings of the 15th ACM SIGPLAN symposium on Principles and practice of parallel programming (PPoPP 10), pp. 1 27 T36 (2010)14. D. Goddeke, R. Strzodka. Cyclic reduction tridiagonal solvers on gpus applied to mixed precision multigrid. IEEE Trans. Parallel Distrib. Syst. 2

23、2, 22 W2 (2011)15. H.S. Stone. Parallel tridiagonal equation solvers. ACM Trans. Math. Softw. 1, 289 4307 (1975)16. H. S. Stone. An efficient parallel algorithm for the solution of a tridiagonal linear system of equations. J. ACM 20, 2738 (1973)17. H. H.Wang. A parallel method for tridiagonal equati

24、ons. ACM Trans. Math. Softw. 7, 170 183 (1981)18. Michielse P, van der Vorst. Data transport in Wangs partition method. Parallel Computing, 7(l): 87 -95 (1988)19. 遲利華，李曉梅.三對(duì)角線性方程組的分布式并行算法.計(jì)算機(jī)研究與發(fā) 展，35(11): 1004-1007 (1998)20. Bondeli S. Divide and conquer: a parallel algorithm for the solution of a

25、tridiagonal linear system of equations. Parallel Comp., 17: 41934 (1991)21. Muller S M Scheerer D. A method to parallelize tridiagonal solvers. Parallel Comput., 17(2&3): 181 -188 (1991)22. Rodrigue G H, Madsen N K, Karush J I. Odd-even reduction for banded linear equations. J.ACM, 26(l): 72-81

26、(1979)23. Meier U. A parallel partition method for solving banded systems of linear equations, Parallel Comput., 2(l): 3343(1985)24. Kapur R N, Browne J C. Techniques for solving block tridiagonal systems on reconfigurable array computers. SIAM J.Sci.Stat. Comput., 5(3): 701719 (1984) 25. RuggieroV,

27、 Galligani E. A parallel algorithm for solving block tridiagonal linear systems, Computers Math. Applic. 24(4): 15-21 (1992)26. 李曉梅，蔣增榮.并行算法.長(zhǎng)沙：湖南科學(xué)技術(shù)出版社,199227. Chung K L, Tsai Y H, Yan W M. A parallel solver for circulant block-tridiagonal systems. Computers Math. Applic, 29(l): 109-113 (1995)28.

28、Evans D J .On the solution of certain Toeplitz tridiagonal linear systems. SIAM J. Numer. Anai., 17(5): 675 -680 (1980)29. Buckley A. On the solution of certain skew symmetric linear systems. SIAM J. Numer. Anal., 14: 566七70 (1977)30. 趙自春，李曉梅.Toeplitz循環(huán)三對(duì)角方程組的并行算法.國(guó)防系統(tǒng)分析 與軟件，2: 50右5 (1989)31. 趙自春，李曉

29、梅.Toeplitz三對(duì)角方程組的兩種并行解法.“第一屆全國(guó)科 學(xué)計(jì)算并行算法論文集”，會(huì)議籌備組編，6: 125 T30 (1989)32. Evans D J. The strides reduction algorithms for solving tridiagonal linear systems. International Journal of Computer Mathematics, 41(3/4): 237250 (1992)33. Evans D J, Yousif W S. The parallel strides reduction algorithms for so

30、lving tridiagonal linear systems. Neural, Parallel and Scientific Computations, 1(l): 29 42 (1993)34. 成禮智，蔣增榮.帶狀(塊)Toeplitz方程組的快速并行算法.數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用，15(1):4451 (1994)35. Buzbee B L, Golub G H, Nielson C W. On direct methods for solving Poisson's equations. SIAM J. Numer. Anal., 7(4): 62755 (1970)36. Buneman O. A compact non-iterative Poisson solver: Report 294 R. Stanford, CA: Stanford University Institute for Plasma Research, 196937. Sweet R A. A cyclic reduction algorithm for solving block tridiagonal systems of arbi