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

邢 茁,呂全義

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

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

邢 茁,呂全義

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

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

線性對(duì)稱變換方程;SYMMLQ算法;并行計(jì)算

0 引 言

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

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

1 SYMMLQ算法的推廣

1.1 求解線性對(duì)稱變換方程的SYMMLQ算法

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

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

T(X)=F,

(1)

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

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

考慮方程(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)

因?yàn)?/p>

(3)

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

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

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

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

(2) 置k:=2,計(jì)算

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

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

1.2 SYMMLQ算法的收斂性分析

證明 由于

因?yàn)镼(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 設(shè)方程(1)相容,則對(duì)任意初始元素X(0)∈V,算法1至多需n步有R(n)=0,即至多需n步算法收斂(n為歐式空間的維數(shù)).

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

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

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

AX+XB=F.

(4)

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

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

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

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

(2) 計(jì)算過程:

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

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

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

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

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

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

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

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

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

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

3.1 算例

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

例1 Poisson方程第一邊值問題

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

當(dāng)n=1 200時(shí),計(jì)算結(jié)果見表1.

表 1 例1計(jì)算結(jié)果

例2 在擋土墻的彈性力學(xué)分析中,由于其結(jié)構(gòu)的形狀和受力、約束情況具有一定的特點(diǎn),所以將應(yīng)力分量之間的關(guān)系近似表示為平衡偏微分方程，即

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

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

計(jì)算結(jié)果見表2.

表 2 例2計(jì)算結(jié)果

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

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

表 3 例3計(jì)算結(jié)果

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

3.2 結(jié)果分析

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

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

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

致謝:感謝西北工業(yè)大學(xué)高性能計(jì)算研究與發(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.矩陣計(jì)算[M].北京:人民郵電出版社,2011.

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

[10] 吳建平,王正華,李曉梅.稀疏線性方程組的高效求解與并行算法[M].長(zhǎng)沙:湖南科學(xué)技術(shù)出版社,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] 蔣爾雄.矩陣計(jì)算[M].北京:科學(xué)出版社,2008.

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

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

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

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

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

[15] 曹方穎,呂全義.預(yù)處理變形共軛梯度法并行求解矩陣的Moore-Penrose 逆[J].紡織高校基礎(chǔ)科學(xué)學(xué)報(bào),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].北京:科學(xué)出版社,2007.

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

[18] 李建江.并行計(jì)算機(jī)及編程基礎(chǔ)[M].北京:清華大學(xué)出版社,2011.

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

編輯、校對(duì)：師 瑯

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

陜西省自然科學(xué)基金資助項(xiàng)目(2009JM1008)

呂全義(1963—),女,遼寧省沈陽市人,西北工業(yè)大學(xué)副教授，研究方向?yàn)樾畔⑻幚碇械目焖倥c并行算法.

E-mail:luquan@nwpu.edu.cn

邢茁,呂全義.一種求解線性對(duì)稱變換方程的并行算法[J].紡織高?；A(chǔ)科學(xué)學(xué)報(bào),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