邢 茁,呂全義

(西北工業(yè)大學 應用數(shù)學系,陜西 西安 710129)

一種求解線性對稱變換方程的并行算法

邢 茁,呂全義

(西北工業(yè)大學 應用數(shù)學系,陜西 西安 710129)

研究求解線性對稱變換方程的SYMMLQ并行算法.將求解線性方程組的SYMMLQ算法推廣應用到求解線性對稱變換方程,將并行過程中的兩次全歸約減少到一次,并對該算法進行改進,以提高并行性,減少計算時間.利用改進后的SYMMLQ算法在并行機上對Poisson方程與橢圓偏微分方程進行效果測試,并與未改進的SYMMLQ算法進行比較和分析.結果表明，改進的SYMMLQ算法的并行效率明顯優(yōu)于未改進的SYMMLQ算法.

線性對稱變換方程;SYMMLQ算法;并行計算

0 引 言

線性變換方程包含最基本的線性方程組、線性矩陣方程,如Lyapunov矩陣方程,Sylvester矩陣方程等.這些方程(組)的求解問題廣泛來源于控制理論[1]、圖像恢復[2]、遷移理論中Riccati微分方程的數(shù)值解[3]、模型降階[4-6]、線性系統(tǒng)的穩(wěn)定性分析[7]等許多科學領域.

目前,求解大型線性變換方程通常采用迭代法,如ADI方法[8](Alternating Direction Implicit)及Krylov子空間法[9-13]等.ADI方法的不便之處在于確定最優(yōu)參數(shù);共軛梯度法[14](conjugate gradient method,簡稱CG法) 具有存儲量小、適合并行計算等優(yōu)點,但它僅適用于系數(shù)矩陣為對稱正定的情形;變形共軛梯度法[15](modified conjugate gradient method,簡稱MCG法)的收斂速度與系數(shù)矩陣的條件數(shù)密切相關,條件數(shù)越小,收斂速度越快,因此MCG算法不適用于條件數(shù)較大的情形; SYMMLQ算法[16]適合求解系數(shù)矩陣為對稱這一情形,其易于實現(xiàn)并行計算且收斂速度快.但是該算法的缺點是在并行實現(xiàn)過程中,需要兩次全歸約,由此引起的大規(guī)模全局通訊使得算法的可擴放性比較差[17-18].

1 SYMMLQ算法的推廣

1.1 求解線性對稱變換方程的SYMMLQ算法

文獻[16]給出了求解對稱情形下線性方程組的SYMMLQ算法,文中將該算法推廣到求解線性對稱變換方程,以拓寬其使用范圍.

設線性變換T是n維歐式空間V上的對稱變換,即?X,Y∈V,有[T(X),Y]=[X,T(Y)].線性對稱變換方程為

T(X)=F,

(1)

其中X∈V是未知的元素,F∈V是已知元素.

根據(jù)文獻[12]的Lanczos過程,推廣到線性對稱變換方程，可得到如下關系式

考慮方程(1)的求解,令

X(k)=X(0)+(Q(1),Q(2),…,Q(k))f(k),

其中f(k)是k×1列向量.取f(k)使得R(k)=F-T(X(k))正交于Q(1),Q(2),…,Q(k),即

[Q(i),R(k)]=0(i=1,…,k).

(2)

因為

(3)

將式(3)左右兩端分別與Q(i)(i=1,…,k)做內積,由式(2)和(3),可得

‖R(0)‖F(xiàn)e1-Tkf(k)=0.

因為對稱三對角矩陣Tk可能奇異,因而對它采用LQ分解,該過程類似于文獻[12].于是,得到求解線性對稱變換方程(1)的SYMMLQ算法(算法1):

(1) 任意給定初始元素X(0)∈V,置k:=1,計算

(2) 置k:=2,計算

(3) 計算ΔX(k-1)=X(k-1)-X(k-2),若‖ΔX(k-1)‖<ε停止計算;否則,計算

(4) 置k:=k+1,轉(3).

1.2 SYMMLQ算法的收斂性分析

證明 由于

因為Q(k+1)正交于Q(1),Q(2),…,Q(k),所以

定理2 殘量‖R(k)‖滿足[R(k),R(j)]=0,其中j=1,2,…,k-1.

證明 由上述定理1的證明過程可知

[Q(k+1),Q(j+1)]=0(j=1,2,…,k-1),

所以[R(k),R(j)]=0. 證畢.

定理3 設方程(1)相容,則對任意初始元素X(0)∈V,算法1至多需n步有R(n)=0,即至多需n步算法收斂(n為歐式空間的維數(shù)).

證明 n維歐式空間V中兩兩正交的非零元素至多為n個,又由定理2的結論可知對任意初始元素X(0)∈V所得到的殘量R(0)與R(1),…,R(n-1),R(n)相互正交,因此算法1至多需n步有R(n)=0,即結論成立. 證畢.

2 SYMMLQ算法的并行實現(xiàn)

針對線性空間V=Rn×n上的對稱變換T(X)=AX+XB,即方程(1)形式如下:

AX+XB=F.

(4)

其中A,B,F∈Rn×n,并且A=AT,B=BT.下面給出求解方程(4)的SYMMLQ算法的具體并行實現(xiàn)過程.

設處理機臺數(shù)為p,p整除n,即n=pl,pi(i=1,2,…,p)表示第i臺處理機.記

注 算法涉及到的所有矩陣乘法皆采用行行劃分矩陣乘法的并行算法[17-18],文中所涉及的矩陣乘法并行計算不再描述.

2.1 未改進的SYMMLQ算法的并行實現(xiàn)

(2) 計算過程:

(3) 循環(huán)過程:

步驟4 在pi中,計算變量ξk=-(εkξk-2+δkξk-1)/rk以及

步驟5 置k:=k+1,返回步驟1.

2.2 改進的SYMMLQ算法的并行實現(xiàn)

在SYMMLQ算法的并行實現(xiàn)過程中,每迭代一次需要計算

ak=[T(Q(k)),Q(k)],bk=‖T(Q(k))-akQ(k)-bk-1Q(k-1)‖，

顯然需要2次全歸約.為了保證改進的SYMMLQ算法的結構穩(wěn)定性,且基本不影響原算法的計算精度,在經過多次嘗試后,對計算步驟進行重新安排,將上述計算步驟調整為如下形式:

這樣,ak,ek的計算只需要1次全歸約,從而達到減少通訊的目的.

調整步驟后,得到了改進的SYMMLQ算法.雖然這個過程需增加一些計算量,但相對很小.同時,改進的算法結構保持不變,所以編寫程序代碼時只需調整原算法相對應的循環(huán)過程步驟2即可,其余步驟保持不變.具體實現(xiàn)過程如下:

3 數(shù)值算例與結果分析

3.1 算例

在HPZ80工作站集群并行機上求解3個算例,分別采用SYMMLQ算法與本文提出的改進的SYMMLQ算法求解,并對計算結果進行比較.在迭代計算過程中,所有的初始矩陣X(0)為零矩陣,終止條件ε=10-6.

例1 Poisson方程第一邊值問題

其中,Γ為單位正方形區(qū)域的邊界,取步長h=1/1201.采用中心差分格式,可將本題轉化為求解矩陣方程AX+XB=F的問題.設矩陣A=(aij)n×n,B=A,F=(fij)n×n,

當n=1 200時,計算結果見表1.

表 1 例1計算結果

例2 在擋土墻的彈性力學分析中,由于其結構的形狀和受力、約束情況具有一定的特點,所以將應力分量之間的關系近似表示為平衡偏微分方程，即

其中,假設正應力u1，u2與切應力u3平衡,即u1=u2=u3=u,并且f(x,y)=x+y,應力邊界條件為u|Γ=0,Γ為受力邊界.于是,平衡偏微分方程轉化為

取步長h=1/1 001,采用中心差分格式以及一些近似處理,可將本題轉化為求解Sylvester矩陣方程AX+XB=F的問題.設矩陣A=(aij)n×n,B=(bij)n×n,F=(fij)n×n,

計算結果見表2.

表 2 例2計算結果

例3 設矩陣A=(aij)n×n,B=(bij)n×n,其中

求解Sylvester矩陣方程AX+XB=F.在此取n=1 500,F為任意給定的矩陣.計算結果見表3.

表 3 例3計算結果

表1～3中p表示處理機臺數(shù),T/s表示時間,K表示迭代次數(shù),E表示相對并行效率,E1表示絕對并行效率,Δ表示誤差‖R(k)‖.由于一臺處理機的計算時間較長,因此采用多臺處理機的計算時間,并且使用多臺處理機中最小的計算時間作為基準來計算加速比和并行效率.

3.2 結果分析

(1) 改進的SYMMLQ算法比未改進的算法減少了并行計算時間,這是由于改進后的算法在循環(huán)迭代中減少了全歸約的次數(shù),從而縮短了全局通訊所引起的時間消耗.

(2) 改進的SYMMLQ算法的迭代次數(shù)與原SYMMLQ算法的迭代次數(shù)基本相同,說明改進的算法沒有破壞原有SYMMLQ算法的結構與數(shù)值穩(wěn)定性,而且由表中誤差數(shù)據(jù)可以看出本文改進的算法對計算的精度影響不大.

(3) 表3數(shù)據(jù)表明當計算規(guī)模一定時,處理機增加到30臺以上,并行效率會急劇下降,這是由于整體內積的計算與全收集而引起的大規(guī)模通訊使得消耗增加,進一步說明本文算法與原算法的可擴放性比較差,因此為保證較高的并行效率應采用適當數(shù)量的處理機.

致謝:感謝西北工業(yè)大學高性能計算研究與發(fā)展中心的大力支持!

[1] DATTA B.Numerical methods for linear control systems[M].Berlin:Elsevier Academic Press,2004.

[2] CALVETTI D,REICHEL L.Application of ADI iterative methods to the restoration of noisy images[J].Siam J Matrix Anal Appl,1996,17(1):165-186.

[3] ENRIGHT W H.Improving the efficiency of matrix operations in the numerical solution of stiff ordinary differential equations[J].ACM Trans Math Software,1978,4(2):127-136.

[4] ANTOULAS A C.Approximation of large-scale dynamical systems (Advances in design and control)[M].Philadelphia:Society of Industnal and Applied Mathematics,2005.

[5] BAUR U,BENNER P.Cross-gramian based model reduction for data-sparse systems[J].Electron Trans Numer Anal,2008,31(1):256-270.

[6] SORENSEN D,ANTOULAS A.The sylvester equation and approximate balanced reduction[J].Linear Algebra Appl,2002,351-352(2): 671-700.

[8] PEACEMAN D W,RACHFORD H H.The numerical solution of parabolic and elliptic differential equations[J].J Soc Ind Appl Math,1995,3(1): 28-41.

[9] GOLUB G H,LOAN C F.矩陣計算[M].北京:人民郵電出版社,2011.

GOLUB G H,LOAN C F.Matrix calculations[M].Beijing: Posts and Telecom Press,2011.

[10] 吳建平,王正華,李曉梅.稀疏線性方程組的高效求解與并行算法[M].長沙:湖南科學技術出版社,2004:269-273.

WU Jianping,WANG Zhenghua,LI Xiaomei.The efficient solution and parallel algorithm of spare linear equations[M].Changsha:Science and Technology Press of Hunan,2004:269-273.

[11] XIE Yajun,HUANG Na,MA Changfeng.Iterative method to solve the generalized coupled Sylvester-transpose linear matrix equations over reflexive or anti-reflexive matrix[J].Computers and Mathematics with Applications,2014,67(11):2071-2084.

[12] 蔣爾雄.矩陣計算[M].北京:科學出版社,2008.

JIANG Erxiong.Matrix calculations[M].Beijing:Science Press,2008.

[13] 李大明.數(shù)值線性代數(shù)[M].北京:清華大學出版社,2010.

LI Daming.Numerical linear algebra[M].Beijing:Tsinghua University Press,2008.

[14] 張凱院,徐仲.數(shù)值代數(shù)[M].2版.北京:科學出版社,2010.

ZHANG Kaiyuan,XU Zhong.Numerical algebra[M].2ed.Beijing:Scince Press,2010.

[15] 曹方穎,呂全義.預處理變形共軛梯度法并行求解矩陣的Moore-Penrose 逆[J].紡織高?；A科學學報,2013,26(1):137-142.

CAO Fangyin,LYU Quanyi.A parallel preconditioned modified conjugate gradient method of Moore-Penrose inverse of matrices[J].Basic Sciences Journal of Textile University,2013,26(1):137-142.

[16] PAIGN C C,SAUNDERS M A.Solution of sparse indefinite system of linear equations[J].Siam J Number Anal,1975,12(4):617-629.

[17] 李曉梅,吳建平.數(shù)值并行算法與軟件[M].北京:科學出版社,2007.

LI Xiaomei,WU Jianping.Numerical parallel algorithm and software[M].Beijing:Science Press,2007.

[18] 李建江.并行計算機及編程基礎[M].北京:清華大學出版社,2011.

LI Jianjiang.Parallel computers and programming base[M].Beijing:Tsinghua University Press,2011.

編輯、校對：師 瑯

A parallel algorithm for solving equation of symmetrical linear transformation

XINGZhuo,LYUQuanyi

(Department of Applied Mathematics, Northwestern Polytechnical University,Xi′an 710129, China)

A parallel algorithm with SYMMLQ method for solving the equation of symmetrical linear transformation was studied.The SYMMLQ method for solving linear equations is extended to solve the equation of symmetrical linear transformation.The degree of reduce operator is decreased from twice to once in the parallel process so as to improve the parallelism of the algorithm and thus reduce the computing time.The Poisson equation and elliptic partial differential equation were tested with the proposed algorithm and the original one, and the results were compared and analyzed.It is shown that the proposed SYMMLQ algorithm is superior to the original SYMMLQ algorithm.

the equation of symmetrical linear transformation;SYMMLQ algorithm;parallel computation

1006-8341(2016)04-0508-07

10.13338/j.issn.1006-8341.2016.04.016

2016-06-11

陜西省自然科學基金資助項目(2009JM1008)

呂全義(1963—),女,遼寧省沈陽市人,西北工業(yè)大學副教授，研究方向為信息處理中的快速與并行算法.

E-mail:luquan@nwpu.edu.cn

邢茁,呂全義.一種求解線性對稱變換方程的并行算法[J].紡織高校基礎科學學報,2016,29(4)：508-514.

XING Zhuo,LYU Quanyi.A parallel algorithm for solving equation of symmetrical linear transformation[J].Basic Sciences Journal of Textile Universities,2016,29(4):508-514.

O 246

A