并行計(jì)算-矩陣特征值計(jì)算--

1、9矩陣特征值計(jì)算在實(shí)際的工程計(jì)算中，經(jīng)常會(huì)遇到求n階方陣A的特征值(Eigenvalue)與特征向量(Eigenvector)的問(wèn)題。對(duì)于一個(gè)方陣A,如果數(shù)值入使方程組Ax=入x即(A-Xn)x=0有非零解向量(SolutionVector)x,則稱(chēng)訥方陣A的特征值，而非零向量x為特征值所對(duì)應(yīng)的特征向量，其中In為n階單位矩陣。由于根據(jù)定義直接求矩陣特征值的過(guò)程比較復(fù)雜，因此在實(shí)際計(jì)算中，往往采取一些數(shù)值方法。本章主要介紹求一般方陣絕對(duì)值最大的特征值的乘哥(Power)法、求對(duì)稱(chēng)方陣特征值的雅可比法和單側(cè)旋轉(zhuǎn)(One-sideRotation)法以及求一般矩陣全部特征值的QR方法及一些相關(guān)的并

2、行算法。1.1 求解矩陣最大特征值的乘哥法1.1.1 乘幕法及其用行算法在許多實(shí)際問(wèn)題中，只需要計(jì)算絕對(duì)值最大的特征值，而并不需要求矩陣的全部特征值。乘哥法是一種求矩陣絕對(duì)值最大的特征值的方法。記實(shí)方陣A的n個(gè)特征值為%i=(1,2,n),且滿(mǎn)足：1213-n特征值入對(duì)應(yīng)的特征向量為xi。乘哥法的做法是：取n維非零向量V0作為初始向量；對(duì)于k=1,2,做如下迭代：直至uk1Uk.,|Uk=AVk-1Vk=Uk/UkJ00卜II的止，這時(shí)Vk+i就是A的絕對(duì)值最大的特征值入i所對(duì)應(yīng)的特征向8量Xi。若Vk-1與Vk的各個(gè)分量同號(hào)且成比例，則入產(chǎn)IUkI若Vk-1與Vk的各個(gè)分量異號(hào)且成比例,則入

3、產(chǎn)-IUkIooo若各取一次乘法和加法運(yùn)算時(shí)間、一次除法運(yùn)算時(shí)間、一次比較運(yùn)算時(shí)間為一個(gè)單位時(shí)間，則因?yàn)橐惠営?jì)算要做一次矩陣向量相乘、一次求最大元操作和一次規(guī)格化操作，所以下述乘哥法串行算法21.1的一輪計(jì)算時(shí)間為n2+2n=O(n2)。算法21.1單處理器上乘哥法求解矩陣最大特征值的算法輸入：系數(shù)矩陣An。，初始向量VnX1,輸出：最大的特征值maxBeginwhile(diff)edo(1)fori=1tondo(1.(1) m=0(1.(2) forj=1tondoSUm=SUm+ai,j*xjendfor(1.(3) i=sumendfor(2)max=b1(3) fori=2tond

4、oif(bimax)thenmax=biendifendfor(4) fori=1tondoxi=bi/maxendfor(5)diff=max-oldmax,oldmax=maxendwhileEnd1.1.2乘幕法并行算法乘哥法求矩陣特征值由反復(fù)進(jìn)行矩陣向量相乘來(lái)實(shí)現(xiàn)，因而可以采用矩陣向量相乘的數(shù)據(jù)劃分方法。設(shè)處理器個(gè)數(shù)為p,對(duì)矩陣A按行劃分為p塊，每塊含有連續(xù)的m行向量，這里m=n/p,編號(hào)為i的處理器含有A的第im至第(i+1)m-1行數(shù)據(jù)，(i=0,1,p-1),初始向量v被廣播給所有處理器。各處理器并行地對(duì)存于局部存儲(chǔ)器中a的行塊和向量v做乘積操作，并求其m維結(jié)果向量中的最大值lo

5、calmax,然后通過(guò)歸約操作的求最大值運(yùn)算得到整個(gè)n維結(jié)果向量中的最大值max并廣播給所有處理器，各處理器利用max進(jìn)行規(guī)格化操作。最后通過(guò)擴(kuò)展收集操作將各處理器中的m維結(jié)果向量按處理器編號(hào)連接起來(lái)并廣播給所有處理器，以進(jìn)行下一次迭代計(jì)算，直至收斂。具體算法框架描述如下：算法21.2乘哥法求解矩陣最大特征值的并行算法輸入：系數(shù)矩陣An。，初始向量Vnx1,e輸出：最大的特征值maxBegin對(duì)所有處理器my_rank(my_rank=0,，p-1)同時(shí)執(zhí)行如下的算法：while(diff)sdo/*diff為特征向量的各個(gè)分量新舊值之差的最大值*/(1)fori=0tom-1do/*對(duì)所存的

6、行計(jì)算特征向量的相應(yīng)分量*/(1.(1) m=0(1.(2) forj=0ton-1dosum=sum+ai,j*xjendfor(1.(3) i=sumendfor(2)localmax=b0/*對(duì)所計(jì)算的特征向量的相應(yīng)分量求新舊值之差的最大值localmax*/(3)fori=1tom-1doif(bi：localmax)thenlocalmax=biendifendfor(4)用Allreduce操作求出所有處理器中l(wèi)ocalmax值的最大值max并廣播到所有處理器中(5)fori=0tom-1do/*對(duì)所計(jì)算的特征向量歸一化*/bi=bi/maxendfor(6)用Allgather操

7、作將各個(gè)處理器中計(jì)算出的特征向量的分量的新值組合并廣播到所有處理器中(7)diff=max-oldmax,oldmax=maxendwhileEnd若各取一次乘法和加法運(yùn)算時(shí)間、一次除法運(yùn)算時(shí)間、一次比較運(yùn)算時(shí)間為一個(gè)單位時(shí)間，一輪迭代的計(jì)算時(shí)間為mn+2m;一輪迭代中，各處理器做一次歸約操作，通信量為1,一次擴(kuò)展收集操作，通信量為m,則通信時(shí)間為4ts(pp-1)十(m+1)tw(p-1)o因此乘哥法的一輪并行計(jì)算時(shí)間為T(mén)p=mn+2m+4ts/-p-1)+(m+1)tw(p-1)。(、MPI源程序請(qǐng)參見(jiàn)所附光盤(pán)。1.2求對(duì)稱(chēng)矩陣特征值的雅可比法1.2.1 雅可比法求對(duì)稱(chēng)矩陣特征值的串行算法

8、若矩陣A=aj是n階實(shí)對(duì)稱(chēng)矩陣，則存在一個(gè)正交矩陣U,使得UTAU=D其中D是一個(gè)對(duì)角矩陣，即入o0D=0屹0I0?九n0這里入(i=1,2,n)是A的特征值，U的第i列是與入對(duì)應(yīng)的特征向量。雅可比算法主要是通過(guò)正交相似變換將一個(gè)實(shí)對(duì)稱(chēng)矩陣對(duì)角化，從而求出該矩陣的全部特征值和對(duì)應(yīng)的特征向量。因此可以用一系列的初等正交變換逐步消去A的非對(duì)角線元素，從而使矩陣A對(duì)角化。設(shè)初等正交矩陣為R(p,q,其(p,p)(q,q)位置的數(shù)據(jù)是cos。，(p,q)(q,p)位置的數(shù)據(jù)分別是-sin酥口sin&pwq),其它位置的數(shù)據(jù)和一個(gè)同階數(shù)的單位矩陣相同。顯然可以得到：R(p,q,)TR(p,q,0)=In

9、不妨記8=R(p,q,9)TAR(p,q,0),如果將右端展開(kāi)，則可知矩陣B的元素與矩陣A的元素之間有如下關(guān)系：bpp=appcos2卅aqqsin20apqsin20bpq=(aqq-app)sin20)/2+apqcos20bpj=apjcos卅aqjsin0bip=aipcos卅aiqsin022八bqq一appsin9+aqqcos9-apqsin20bqp=bpqbqj=-apjsin0+aqjcos0biq=-aipsin0+aiqcos0bij=aijR為正交矩陣，所以B亦為對(duì)稱(chēng)矩陣。若其中1wi,jwn且i,jwp,q。因?yàn)锳為對(duì)稱(chēng)矩陣,要求上I陣B的元素bpq=0,則只需令(

10、aqq-app)sin20)/2+apqcos29=0,即:pqtg2u=(aqq-app)2sin20,sin0和cos0就可以了。在實(shí)際應(yīng)用時(shí)，考慮到并不需要解出&而只需要求出如果限制。值在-兀/220)do(1.(1) max=a1,2(1.(2) fori=1tondoforj=i+1tondoif(ai，j)max)thenmax=ai，j,p=i,q=jendifendforendfor(1.(3) mpute:m=-ap,q,n=(aq,q-ap,p)/2,w=sgn(n)*m/sqrt(m*m+n*n),si2=w,si1=w/sqrt(2(1+sqrt(1-w*w),co1=

11、sqrt(1-si1*si1),bp,p=ap,p*co1*co1+aq,q*si1*si1+ap,q*si2,bq,q=ap,p*si1*si1+aq,q*co1*co1-ap,q*si2,bq,p=0,bp,q=0(1.(4) forj=1tondoif(jwp)and(jwq)then(i)bp,j=ap,j*co1+aq,j*si1(ii)bq,j=-ap，j*si1+aq,j*co1endifendfor(1.(5) fori=1tondoEndif(iwp)and(iwq)then(i)bi,p=ai,p*co1+ai,q*si1(11) bi,q=-ai,p*si1+ai,q*c

12、o1endifendfor(1.(6) fori=1tondoforj=1tondoai,j=bi,jendforendforendwhile(2)fori=1tondoEigenvaluei=ai,iendfor1.2.2雅可比法求對(duì)稱(chēng)矩陣特征值的并行算法串行雅可比算法逐次尋找非主對(duì)角元絕對(duì)值最大的元素的方法并不適合于并行計(jì)算。此，在并行處理中，我們每次并不尋找絕對(duì)值最大的非主對(duì)角元消去，而是按一定的順序?qū)⒅械乃猩先窃厝肯ヒ槐椋@樣的過(guò)程稱(chēng)為一輪。由于對(duì)稱(chēng)性，在一輪中，A的三角元素也同時(shí)被消去一遍。經(jīng)過(guò)若干輪，可使A的非主對(duì)角元的絕對(duì)值減少，收斂于個(gè)對(duì)角方陣。具體算法如下：設(shè)處理器

13、個(gè)數(shù)為p對(duì)矩陣A按行劃分為2P塊每塊含有連續(xù)的m/2行向量，記m=n/p,些行塊依次記為Ao,Ai,A2P，并將A2i與A2i+1存放與標(biāo)號(hào)為i的處理器中。每輪計(jì)算開(kāi)始，各處理器首先對(duì)其局部存儲(chǔ)器中所存的兩個(gè)行塊的所有行兩兩配對(duì)進(jìn)行旋轉(zhuǎn)變換，消去相應(yīng)的非對(duì)角線元素。然后按圖21.1所示規(guī)律將數(shù)據(jù)塊在不同處理器之間傳送，以消去其它非主對(duì)角元素。圖1.1p=4時(shí)的雅可比算法求對(duì)稱(chēng)矩陣特征值的數(shù)據(jù)交換圖這里長(zhǎng)方形框中兩個(gè)方格內(nèi)的整數(shù)被看成是所移動(dòng)的行塊的編號(hào)。在要構(gòu)成新的行塊配對(duì)時(shí)，只要將每個(gè)處理器所對(duì)應(yīng)的行塊按箭頭方向移至相鄰的處理器即可，這樣的傳送可以在行塊之間實(shí)現(xiàn)完全配對(duì)。當(dāng)編號(hào)為i和j的兩個(gè)

14、行塊被送至同一處理器時(shí)，令編號(hào)為i的行塊中的每行元素和編號(hào)為j的行塊中的每行元素配對(duì)，以消去相應(yīng)的非主對(duì)角元素，這樣在所有的行塊都兩兩配對(duì)之后，可以將所有的非主對(duì)角元素消去一遍，從而完成一輪計(jì)算。由圖1.1可以看出，在一輪計(jì)算中，處理器之間要2p-1次交換行塊。為保證不同行塊配對(duì)時(shí)可以將原矩陣A的非主對(duì)角元素消去，引入變量b記錄每個(gè)處理器中的行塊數(shù)據(jù)在原矩陣A中的實(shí)際行號(hào)。在行塊交換時(shí)，變量b也跟著相應(yīng)的行塊在各處理器之間傳送。處理器之間的數(shù)據(jù)塊交換存在如下規(guī)律：0號(hào)處理器前一個(gè)行塊(簡(jiǎn)稱(chēng)前數(shù)據(jù)塊，后一個(gè)行塊簡(jiǎn)稱(chēng)后數(shù)據(jù)塊)始終保持不變，將后數(shù)據(jù)塊發(fā)送給1號(hào)處理器，作為1號(hào)處理器的前數(shù)據(jù)塊。同時(shí)

15、接收1號(hào)處理器發(fā)送的后數(shù)據(jù)塊作為自己的后數(shù)據(jù)塊。p-1號(hào)處理器首先將后數(shù)據(jù)塊發(fā)送給編號(hào)為p-2的處理器，作為p-2號(hào)處理器的后數(shù)據(jù)塊。然后將前數(shù)據(jù)塊移至后數(shù)據(jù)塊的位置上，最后接收p-2號(hào)處理器發(fā)送的前數(shù)據(jù)塊作為自己的前數(shù)據(jù)塊。編號(hào)為my_rank的其余處理器將前數(shù)據(jù)塊發(fā)送給編號(hào)為my_rank+1的處理器，作為my_rank+1號(hào)處理器的前數(shù)據(jù)塊。將后數(shù)據(jù)塊發(fā)送給編號(hào)為my_rank-1的處理器，作為my_rank-1號(hào)處理器的后數(shù)據(jù)塊。為了避免在通信過(guò)程中發(fā)生死鎖，奇數(shù)號(hào)處理器和偶數(shù)號(hào)處理器的收發(fā)順序被錯(cuò)開(kāi)。使偶數(shù)號(hào)處理器先發(fā)送后接收；而奇數(shù)號(hào)處理器先將數(shù)據(jù)保存在緩沖區(qū)中，然后接收偶數(shù)號(hào)處理

16、器發(fā)送的數(shù)據(jù)，最后再將緩沖區(qū)中的數(shù)據(jù)發(fā)送給偶數(shù)號(hào)處理器。數(shù)據(jù)塊傳送的具體過(guò)程描述如下：算法21.4雅可比法求對(duì)稱(chēng)矩陣特征值的并行過(guò)程中處理器之間的數(shù)據(jù)塊交換算法輸入：矩陣A的行數(shù)據(jù)塊和向量b的數(shù)據(jù)塊分布于各個(gè)處理器中輸出：在處理器陣列中傳送后的矩陣A的行數(shù)據(jù)塊和向量b的數(shù)據(jù)塊Begin對(duì)所有處理器my_rank(my_rank=0,，p-1)同時(shí)執(zhí)行如下的算法：/*矩陣A和向量b為要傳送的數(shù)據(jù)塊*/(1) if(my-rank=0)then/*0號(hào)處理器*/(1.(1) 數(shù)據(jù)塊發(fā)送給1號(hào)處理器(1.(2) 1號(hào)處理器發(fā)送來(lái)的后數(shù)據(jù)塊作為自己的后數(shù)據(jù)塊endif(2) if(my-rank=p-

17、1)and(my-rankmod2豐0then/*處理器p-1且其為奇數(shù)*/(2.1)fori=m/2tom-1do/*將后數(shù)據(jù)塊保存在緩沖區(qū)buffer中*/forj=0ton-1dobufferi-m/2,j=ai,jendforendfor(2.(2) fori=m/2tom-1dobufi-m/2=biendfor(2.(3) fori=0tom/2-1do/*將前數(shù)據(jù)塊移至后數(shù)據(jù)塊的位置上*/forj=0ton-1doai+m/2,j=ai,jendforendfor(2.(4) fori=0tom/2-1dobi+m/2=biendfor(2.(5) p-2號(hào)處理器發(fā)送的數(shù)據(jù)塊作為

18、自己的前數(shù)據(jù)塊(2.(6) uffer中的后數(shù)據(jù)塊發(fā)送給編號(hào)為p-2的處理器endifif(my-rank=p-1)and(my-rankmod2=0)then/*處理器p-1且其為偶數(shù)*/(3.(1) 數(shù)據(jù)塊發(fā)送給編號(hào)為p-2的處理器(3.(2) fori=0tom/2-1do/*將前數(shù)據(jù)塊移至后數(shù)據(jù)塊的位置上*/forj=0ton-1doai+m/2,j=ai,jendforendfor(3.(3) fori=0tom/2-1dobi+m/2=biendfor(3.(4) p-2號(hào)處理器發(fā)送的數(shù)據(jù)塊作為自己的前數(shù)據(jù)塊endif(4) if(my-rankwp-1)and(my-rank豐0

19、then/*其它的處理器*/(4.(1) if(my-rankmod2=0)then/*偶數(shù)號(hào)處理器*/(i)將前數(shù)據(jù)塊發(fā)送給編號(hào)為my_rank+1的處理器(ii)將后數(shù)據(jù)塊發(fā)送給編號(hào)為my_rank-1的處理器(ii)接收編號(hào)為my_rank-1的處理器發(fā)送的數(shù)據(jù)塊作為自己的前數(shù)據(jù)塊(iv)接收編號(hào)為my_rank+1的處理器發(fā)送的數(shù)據(jù)塊作為自己的后數(shù)據(jù)塊else/*奇數(shù)號(hào)處理器*/(v)fori=0tom-1do/*將前后數(shù)據(jù)塊保存在緩沖區(qū)buffer中*/forj=0ton-1dobufferi,j=ai,jendforendfor(vi)fori=0tom-1dobufi=biend

20、for(vii)接收編號(hào)為my_rank-1的處理器發(fā)送的數(shù)據(jù)塊作為自己的前數(shù)據(jù)塊Endendif(viii)接收編號(hào)為my_rank+1的處理器發(fā)送的數(shù)據(jù)塊作為自己的后數(shù)據(jù)塊(ix)將存于buffer中的前數(shù)據(jù)塊發(fā)送給編號(hào)為my_rank+1的處理器(x)將存于buffer中的后數(shù)據(jù)塊發(fā)送給編號(hào)為my_rank-1的處理器endif各處理器并行地對(duì)其局部存儲(chǔ)器中的非主對(duì)角元素aj進(jìn)行消去，首先計(jì)算旋轉(zhuǎn)參數(shù)并對(duì)第i行和第j行兩行元素進(jìn)行旋轉(zhuǎn)行變換。然后通過(guò)擴(kuò)展收集操作將相應(yīng)的旋轉(zhuǎn)參數(shù)及第i列和第j列的列號(hào)按處理器編號(hào)連接起來(lái)并廣播給所有處理器。各處理器在收到這些旋轉(zhuǎn)參數(shù)和列號(hào)之后，按0,1，,

21、p-1的順序依次讀取旋轉(zhuǎn)參數(shù)及列號(hào)并對(duì)其m行中白第i列和第j列元素進(jìn)行旋轉(zhuǎn)列變換。經(jīng)過(guò)一輪計(jì)算的2p-1次數(shù)據(jù)交換之后，原矩陣A的所有非主對(duì)角元素都被消去一次。此時(shí)，各處理器求其局部存儲(chǔ)器中的非主對(duì)角元素的最大元localmax,然后通過(guò)歸約操作的求最大值運(yùn)算求得將整個(gè)n階矩陣非主對(duì)角元素的最大元max,并廣播給所有處理器以決定是否進(jìn)行下一輪迭代。具體算法框架描述如下：算法21.5雅可比法求對(duì)稱(chēng)矩陣特征值的并行算法輸入：矩陣Anxn,輸出：矩陣A的主對(duì)角元素即為A的特征值Begin對(duì)所有處理器my_rank(my_rank=0,，p-1)同時(shí)執(zhí)行如下的算法：(a)fori=0tom-1dobi

22、=myrank*m+i/*b記錄處理器中的行塊數(shù)據(jù)在原矩陣A中的實(shí)際行號(hào)*/endfor(b)while(|max|)edo/*max為A中所有非對(duì)角元最大的絕對(duì)值*/(1)fori=my_rank*mto(my_rank+1)*m-2do/*對(duì)本處理器內(nèi)部所有行兩兩配對(duì)進(jìn)行旋轉(zhuǎn)變換*/forj=i+1to(my_rank+1)*m-1do(1.(1) imodm,t=jmodm/*i,j為進(jìn)行旋轉(zhuǎn)變換行(稱(chēng)為主行)的實(shí)際行號(hào),r,t為它們?cè)趬K內(nèi)的相對(duì)行號(hào)*/(1.(2) (ar,jw0)then/*對(duì)四個(gè)主元素的旋轉(zhuǎn)變換*/(i)Compute:f=-ar,j,g=(at,j-ar,i)/2

23、,h=sgn(g)*f/sqrt(f*f+g*g),sin2=h,sin1=h/sqrt(2*(1+sqrt(1-h*h),cos1=sqrt(1-sin1*sin1),bpp=ar,i*cos1*cos1+at,j*sin1*sin1+ar,j*sin2,bqq=ar,i*sin1*sin1+at,j*cos1*cos1-ar,j*sin2,bpq=0,bqp=0(ii) forv=0ton-1do/*對(duì)兩個(gè)主行其余元素的旋轉(zhuǎn)變換*/if(vwi)and(vwj)thenbrv=ar,v*cos1+at,v*sin1at,v=-ar,v*sin1+at,v*cos1endifendfor(i

24、ii) forv=0ton-1doif(vwi)and(vwj)thenar,v=brvendifendfor(iv) forv=0tom-1do/*對(duì)兩個(gè)主列在本處理器內(nèi)的其余元素的旋轉(zhuǎn)變換*/if(vwr)and(vwt)thenbiv=av,i*cosl+av,j*sin1av,j=-av,i*sin1+av,j*cos1endifendfor(v) forv=0tom-1doif(vwr)and(vwt)thenav,i=bivendifendfor(vi)Compute:ar,i=bpp,at,j=bqq,ar,j=bpq,at,i=bqp,/*用temp1保存本處理器主行的行號(hào)和旋

25、轉(zhuǎn)參數(shù)*/temp10=sin1,temp11=cos1,temp12=(float)i,temp13=(float)jelse(vii)Compute:temp10=0,temp11=0,temp12=0,temp13=0endif(1.(3) 有處理器temp1中的旋轉(zhuǎn)參數(shù)及主行的行號(hào)按處理器編號(hào)連接起來(lái)并廣播給所有處理器，存于temp2中(1.(4) rrent=0(1.(5) forv=1topdo/*根據(jù)temp2中的其它處理器的旋轉(zhuǎn)參數(shù)及主行的行號(hào)對(duì)相關(guān)的列在本處理器的部分進(jìn)行旋轉(zhuǎn)變換*/(i)Compute:s1=temp2(v-1)*4+0,c1=temp2(v-1)*4+1,

26、i1=(int)temp2(v-1)*4+2,j1=(int)temp2(v-1)*4+3(ii)if(s1、c1、i1、j1中有一不為0)thenif(my-rankcurrent)thenforz=0tom-1doziz=az,i1*c1+az,j1*s1az,j1=-az,i1*s1+az,j1*c1endforforz=0tom-1doaz,i1=zizendforendifendif(111) current=current+1endforendforendfor(2) forcounter=1to2p-2do/*進(jìn)彳T2p-2次處理器間的數(shù)據(jù)交換，并對(duì)交換后處理器中所有行兩兩配對(duì)進(jìn)

27、行旋車(chē)t變換*/(2.(1) ta_exchange()/*處理器間的數(shù)據(jù)交換*/(2.(2) fori=0tom/2-1doforj=m/2tom-1do(1) if(ai,bjw0)en/*對(duì)四個(gè)主元素的旋轉(zhuǎn)變換*/Compute:f=-ai,bj,g=(aj,bj-ai,bi)/2,h=sgn(g)*f/sqrt(f*f+g*g),sin2=h,sin1=h/sqrt(2*(1+sqrt(1-h*h),cos1=sqrt(1-sin1*sin1),bpp=ai,bi*cos1*cos1+aj,bj*sin1*sin1+ai,bj*sin2,bqq=ai,bi*sin1*sin1+aj,b

28、j*cos1*cos1-ai,bj*sin2,bpq=0,bqp=0forv=0ton-1do/*對(duì)兩個(gè)主行其余元素的旋轉(zhuǎn)變換*/if(vwbi)and(vwbj)thenbrv=ai,v*cos1+aj,v*sin1aj,v=-ai,v*sin1+aj,v*cos1endifendfor forv=0ton-1doif(vwbi)and(vwbj)thenai,v=brvendifendfor forv=0tom-1do/*對(duì)本處理器內(nèi)兩個(gè)主列的其余元素旋轉(zhuǎn)變換*/if(vwi)and(vwj)thenbiv=av,bi*cos1+av,bj*sin1av,bj=-av,bi*sin1+av

29、,bj*cos1endifendfor forv=0tom-1doif(vwi)and(v刊)thenav,bi=bivendifendfor Compute:ai,bi=bpp,aj,bj=bqq,ai,bj=bpq,aj,bi=bqp/*用tempi保存本處理器主行的行號(hào)和旋轉(zhuǎn)參數(shù)*/temp10=sini,temp11=cosi,temp12=(float)bi,temp13=(float)bjendforelse Compute:temp10=0,temp11=0,temp12=0,temp13=0endif(ii)將所有處理器temp1中的旋轉(zhuǎn)參數(shù)及主行的行號(hào)按處理器編號(hào)連接起來(lái)并廣

30、播給所有處理器，存于temp2中(iii) current=0(iv) forv=1topdo/*根據(jù)temp2中的其它處理器的旋轉(zhuǎn)參數(shù)及主行的行號(hào)對(duì)相關(guān)的列在本處理器的部分進(jìn)行旋轉(zhuǎn)變換*/ Compute:s1=temp2(v-1)*4+0,c1=temp2(v-1)*4+1,i1=(int)temp2(v-1)*4+2,j1=(int)temp2(v-1)*4+3 if(s1、c1、i1、j1中有一不為0)thenif(my-rankcurrent)thenforz=0tom-1doziz=az,i1*c1+az,j1*s1az,j1=-az,i1s1+az,j1*c1endforforz

31、=0tom-1doaz,i1=zizendforendifendif current=current+1endforendforendfor(3)Data-exchange()/*進(jìn)行一輪中的最后一次處理器間的數(shù)據(jù)交換，使數(shù)據(jù)回到原來(lái)的位置*/(4)localmax=max/*localmax為本處理器中非對(duì)角元最大的絕對(duì)值*/(5)fori=0tom-1doforj=0ton-1doif(m*my-rank+i)矜thenEndif(ai,jlocalmax)thenlocalmax=ai,jendifendifendforendfor(6)通過(guò)Allreduce操作求出所有處理器中l(wèi)oca

32、lmax的最大值為新的max值endwhile在上述算法中，每個(gè)處理器在一輪迭代中要處理2p個(gè)行塊對(duì)，由于每一個(gè)行塊對(duì)含有m/2行，因而對(duì)每一個(gè)行塊對(duì)的處理要有(m/2)2=m2/4個(gè)行的配對(duì)，即消去m2/4個(gè)非主對(duì)角元素.每消去一個(gè)非主對(duì)角元素要對(duì)同行n個(gè)元素和同列n個(gè)元素進(jìn)行旋轉(zhuǎn)變換.由于按行劃分?jǐn)?shù)據(jù)塊，對(duì)同行n個(gè)元素進(jìn)行旋轉(zhuǎn)變換的過(guò)程是在本處理器中進(jìn)行的.而對(duì)同列n個(gè)元素進(jìn)行旋轉(zhuǎn)變換的過(guò)程則分布在所有處理器中進(jìn)行.但由于所有處理器同時(shí)在為其它處理器的元素消去過(guò)程進(jìn)行列的旋轉(zhuǎn)變換，對(duì)每個(gè)處理器而言，對(duì)列元素進(jìn)行旋轉(zhuǎn)變換的處理總量仍然是n個(gè)元素。因此，一個(gè)處理器每消去一個(gè)非主對(duì)角元素共要對(duì)2

33、n個(gè)元素進(jìn)行旋轉(zhuǎn)變換。而對(duì)一個(gè)元素進(jìn)行旋轉(zhuǎn)變換需要2次乘法和1次加法，若取一次乘法運(yùn)算時(shí)間或一次加法運(yùn)算時(shí)間為一個(gè)單位時(shí)間，則其需要3個(gè)單位運(yùn)算時(shí)間，即一輪迭代的計(jì)算時(shí)間為T(mén)i=3X2pX2nxm2/4=3n3/p；在每輪迭代中，各個(gè)處理器之間以點(diǎn)對(duì)點(diǎn)的通信方式相互錯(cuò)開(kāi)交換數(shù)據(jù)2p-1次，通信量為mn+m,擴(kuò)展收集操作n(n-1)/(2p)次，通信量為4,另外有歸約操作1次，通信量為1，從而不難得出上述算法求解過(guò)程中的總通信時(shí)間為：T2=4ts+m(n+1)tw(4p-2)+n(n-1)/p+2t3p-1)+2n(n-1)/p+1tw(p-1)(因此雅可比算法求對(duì)矩陣特征值的一輪并行計(jì)算時(shí)間為

34、Tp=T1+T2。我們的大量實(shí)驗(yàn)結(jié)果說(shuō)明，對(duì)于n階方陣，用上述算法進(jìn)行并行計(jì)算，一般需要不超過(guò)O(logn)輪就可以收斂。MPI源程序請(qǐng)參見(jiàn)章末附錄。1.3求對(duì)稱(chēng)矩陣特征值的單側(cè)旋轉(zhuǎn)法1.3.1 單側(cè)旋轉(zhuǎn)法的算法描述求解對(duì)稱(chēng)方陣特征值及特征向量的雅可比算法在每次消去計(jì)算前都要對(duì)非主對(duì)角元素選擇最大元，導(dǎo)致非實(shí)際計(jì)算開(kāi)銷(xiāo)很大。在消去計(jì)算時(shí)，必須對(duì)方陣同時(shí)進(jìn)行行、列旋轉(zhuǎn)變換，這稱(chēng)之為雙側(cè)旋轉(zhuǎn)(Two-sideRotation)變換。在雙側(cè)旋轉(zhuǎn)變換中，方陣的行、列方向都有數(shù)據(jù)相關(guān)關(guān)系，使得整個(gè)計(jì)算中的相關(guān)關(guān)系復(fù)雜，特別是在對(duì)高階矩陣進(jìn)行特征值求解時(shí)，增加了處理器間的通信開(kāi)銷(xiāo)。針對(duì)傳統(tǒng)雅可比算法的缺點(diǎn)

35、，可以將雙側(cè)旋轉(zhuǎn)改為單側(cè)旋轉(zhuǎn)，得出一種求對(duì)稱(chēng)矩陣特征值的快速算法。這一算法對(duì)矩陣僅實(shí)施列變換，使得數(shù)據(jù)相關(guān)關(guān)系僅在同列之間，因此方便數(shù)據(jù)劃分，適合并行計(jì)算，稱(chēng)為單側(cè)旋轉(zhuǎn)法(One-sideRotation)若A為一對(duì)稱(chēng)方陣，入是A的特征值，x是A的特征向量，則有：Ax=及。又A=AT,所以ATAx=Ax=Xx,這說(shuō)明了是ATA的特征值，x是ATA的特征向量，即ATA的特征值是A的特征值的平方，并且它們的特征向量相同。我們使用18.7.1節(jié)中所介紹的Givens旋轉(zhuǎn)變換對(duì)A進(jìn)行一系列的列變換，得到方陣Q使其各列兩兩正交，即AV=Q這里V為表示正交化過(guò)程的變換方陣。因QtQ=(AV)TAV=VTA

36、TAV=diag(B1,B2,，Sn)為n階對(duì)角方陣，可見(jiàn)這里B1,62,，Bn是矩陣ATA的特征值，VI各列是Ata的特征向量。由此可得：6i=&(i=1,2,n)。設(shè)Q的第i列為q,則iqiTq=令腳l|2=e)p=0fori=1tondoforj=i+1tondo(1.1)sum0=0,sum1=0,sum2=0(1.2)fork=1tondosum0=sum0+ak,i*ak,jsum1=sum1+ak,i*ak,isum2=sum2+ak,j*ak,jendfor(1.3)if(sum0e)then(I) aa=2*sum0bb=sum1-sum2rr=sqrt(aa*aa+bb*b

37、b)(II) if(bb0)thenc=sqrt(bb+rr)/(2*rr)s=aa/(2*rr*c)elses=sqrt(rr-bb)/(2*rr)endifc=aa/(2*rr*s)(III) fork=1tondotempk=c*ak,i+s*ak,jak,j=-s*ak,i+c*ak,jendfor(IV) fork=1tondotemp1k=c*ek,i+s*ek,jek,j=-s*ek,i+c*ek,jendfor(v)fork=1tondoendforendforendwhile(2)fori=1tondosu=0ak,i=tempkendfor(vi) fork=1tondoe

38、k,i=temp1kendfor(vii) if(sum0p)thenp=sum0endifendifEndforj=1tondosu=su+aj,i*aj,iendforBj=sqrtsuif(e0,j*a0,j=0)thenc=sqrt(bb+rr)/(2*rr)s=aa/(2*rr*c)elses=sqrt(rr-bb)/(2*rr)c=aa/(2*rr*s)endif forv=0toq-1dotempv=c*av,i+s*av,jav,j=(-s)*av,i+c*av,jendfor forv=0toz-1dotemp1v=c*ev,i+s*ev,jev,j=(-s)*ev,i+c*

39、ev,jendfor forv=0toq-1doav,i=tempvendfor forv=0toz-1doev,i=temp1vendfor k=k+1endifendforendforendwhile0號(hào)處理器從其余各處理器中接收子矩陣a和e,得到經(jīng)過(guò)變換的矩陣A和I:/*0號(hào)處理器負(fù)責(zé)計(jì)算特征值*/(b)fori=1tondoEnd(1)su=0(2) forj=1tondosu=su+Aj,i*Aj,iendfor(3)Bj=sqrt(su)(4)if(I0,j*a0,j)and(countp)thenp=ai,jendifendforEndendforendwhile(2) fori

40、=1tondoEigenvaluei=ai,iendfor顯然，QR分解求矩陣特征值算法的一輪時(shí)間復(fù)雜度為QR分解、矩陣轉(zhuǎn)置和矩陣相乘算法的時(shí)間復(fù)雜度之和。1.4.2 QR方法求一般矩陣全部特征值的并行算法并行QR分解求矩陣特征值的思想就是反復(fù)運(yùn)用并行QR分解和并行矩陣相乘算法進(jìn)行迭代，直到矩陣序列a收斂于一個(gè)上三角矩陣或塊上三角矩陣為止。具體的并行算法描述如下：算法21.9QR方法求一般矩陣全部特征值的并行算法輸入：矩陣An。,單位矩陣Q,輸出：矩陣特征值日genvalueBegin對(duì)所有處理器my_rank(my_rank=0,，p-1)同時(shí)執(zhí)行如下的算法：(a)while(p)and(c

41、ount0)and(my_rankp)thenp=ai,jendifendforendforendwhile(b)fori=1tondoEigenvalue。=ai,iendfor在上述算法的一輪迭代中，實(shí)際上是使用算法18.12、18.1、18.6并行地進(jìn)行矩陣的QR分解、矩陣的轉(zhuǎn)置、矩陣的相乘各一次，因此，一輪迭代的計(jì)算時(shí)間為上述算法的時(shí)間復(fù)雜度之和。MPI源程序請(qǐng)參見(jiàn)所附光盤(pán)。1.5小結(jié)矩陣的特征值與特征向量在工程技術(shù)上應(yīng)用十分廣泛，在諸如系統(tǒng)穩(wěn)定性問(wèn)題和數(shù)字信號(hào)處理的K-L變換中，它都是最基本、常用的算法之一。本章主要討論了矩陣特征值的幾種基本算法，關(guān)于該方面的其它講解與分析讀者可參考文獻(xiàn)1、2。參考文獻(xiàn)1 .陳國(guó)良編著.并行算法的設(shè)計(jì)與分析（修訂版）.高等教育出版社，2002.112 .陳國(guó)良，陳峻編著.VLSI計(jì)算理論與并行算法.中國(guó)科學(xué)技術(shù)大學(xué)出版社，1991附錄求對(duì)稱(chēng)矩陣特征值的雅可比并行算法MPI源程序源程序cjacobi.c#includestdio.h#includestdlib.h#includempi.h#includemath.h#includestring.h#defineE0.00