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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

16、法( recursive doubling)和循環(huán)遞減算法(cyclic reduction)同樣是備受關(guān)注的典型并行方法，已經(jīng)被多種并行計算技術(shù) 實現(xiàn)。塊(循環(huán))三對角線性方程組的求解同樣也是諸多應(yīng)用問題的重要組成部 分。例如，周期的樣條插值就導(dǎo)致對角占優(yōu)的循環(huán)三對角線性方程組的求解，當邊界條件為周期邊界條件時，一些偏微分方程的離散化也可能導(dǎo)致循環(huán)三對角線 性方程組的求解。4參考文獻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. 張寶琳，谷同祥，莫則堯.數(shù)值并行計算原理與方法.北京:國防工業(yè)出版 社，19994.

18、 蔡放，方透，李彬.具有完全并行度的SOR算法.高等學(xué)校計算數(shù)學(xué)學(xué)報， 1:11 T4 (2002)5. 呂全義，葉天麒.系數(shù)矩陣為塊三對角的線性方程組的并行算法.西北工業(yè)大學(xué)學(xué)報,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. 肖曼玉，呂全義.塊三對角線性方程組的一種有效并行算法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. 遲利華，李曉梅.三對角線性方程組的分布式并行算法.計算機研究與發(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. 李曉梅，蔣增榮.并行算法.長沙：湖南科學(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)三對角方程組的并行算法.國防系統(tǒng)分析 與軟件，2: 50右5 (1989)31. 趙自春，李曉

29、梅.Toeplitz三對角方程組的兩種并行解法.“第一屆全國科 學(xué)計算并行算法論文集”，會議籌備組編，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ù)值計算與計算機應(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