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

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

2、值的 QR 方法及 一些相關(guān)的并行算法。1.1 求解矩陣最大特征值的乘冪法1.1.1 乘冪法及其串行算法在許多實(shí)際問(wèn)題中，只需要計(jì)算絕對(duì)值最大的特征值，而并不需要求矩陣的全部特征值。 乘冪法是一種求矩陣絕對(duì)值最大的特征值的方法。記實(shí)方陣A的n個(gè)特征值為i i=(1,2, ,n)， 且滿足：1 2 3 n 特征值i 對(duì)應(yīng)的特征向量為xi 。乘冪法的做法是：取n維非零向量v0 作為初始向量；對(duì)于k=1,2, ,做如下迭代：直至 uk +1 ¥ - ukuk =Avk-1vk = uk /uk < 為止，這時(shí)vk+1 就是A的絕對(duì)值最大的特征值1 所對(duì)應(yīng)的特征向¥量x1 。

3、若vk-1 與vk 的各個(gè)分量同號(hào)且成比例，則1 =uk ；若vk-1 與vk 的各個(gè)分量異號(hào)且成比 例，則1 = -uk 。若各取一次乘法和加法運(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×n ，初始向量v n×1 ，輸出：最大的特征值 maxBeginwhile (diff>) do(1)for i=1 to n do(1.1)s

4、um=0(1.2)for j= 1 to n do sum=sum+ai,j*xj end for(1.3)bi= sumend for (2)max=b1 (3)for i=2 to n doif (bi>max) then max=bi end if end for(4)for i=1 to n doxi =bi/maxend for(5)diff=max-oldmax , oldmax=maxend whileEnd1.1.2 乘冪法并行算法乘冪法求矩陣特征值由反復(fù)進(jìn)行矩陣向量相乘來(lái)實(shí)現(xiàn)，因而可以采用矩陣向量相乘的數(shù) 據(jù)劃分方法。設(shè)處理器個(gè)數(shù)為 p，對(duì)矩陣 A 按行劃分為 p 塊，

5、每塊含有連續(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é)果 向量中的最大值 localmax，然后通過(guò)歸約操作的求最大值運(yùn)算得到整個(gè) n 維結(jié)果向量中的最 大值 max 并廣播給所有處理器，各處理器利用 max 進(jìn)行規(guī)格化操作。最后通過(guò)擴(kuò)展收集操 作將各處理器中的 m 維結(jié)果向量按處理器編號(hào)連接起來(lái)并廣播給所有處理器，以進(jìn)行下一 次迭代計(jì)算，直至收斂。

6、具體算法框架描述如下：算法 21.2乘冪法求解矩陣最大特征值的并行算法輸入：系數(shù)矩陣An×n ，初始向量v n×1 ，輸出：最大的特征值 maxBegin對(duì)所有處理器 my_rank(my_rank=0, p-1)同時(shí)執(zhí)行如下的算法:while (diff>) do/* diff 為特征向量的各個(gè)分量新舊值之差的最大值*/ (1)for i=0 to m-1 do/*對(duì)所存的行計(jì)算特征向量的相應(yīng)分量*/(1.1)sum=0(1.2)for j= 0 to n-1 dosum=sum+ai,j*xjend for(1.3)bi= sumend for(2)localma

7、x=b0/*對(duì)所計(jì)算的特征向量的相應(yīng)分量 求新舊值之差的最大值 localmax */(3)for i=1 to m-1 doif (bi>localmax) then localmax=bi end if end for(4)用 Allreduce 操作求出所有處理器中 localmax 值的最大值 max并廣播到所有處理器中(5)for i=0to m-1 do/*對(duì)所計(jì)算的特征向量歸一化 */bi =bi/maxend for(6)用 Allgather 操作將各個(gè)處理器中計(jì)算出的特征向量的分量的新值組合并廣播到 所有處理器中(7)diff=max-oldmax, oldmax=m

8、axend while End若各取一次乘法和加法運(yùn)算時(shí)間、一次除法運(yùn)算時(shí)間、一次比較運(yùn)算時(shí)間為一個(gè)單位時(shí)間，一輪迭代的計(jì)算時(shí)間為 mn+2m；一輪迭代中，各處理器做一次歸約操作，通信量為 1，一次擴(kuò)展收集操作，通信量為 m，則通信時(shí)間為 4t s (p - 1) + (m + 1)t w ( p - 1) 。因此乘冪法的一輪并行計(jì)算時(shí)間為T p = mn + 2m + 4t s (MPI 源程序請(qǐng)參見(jiàn)所附光盤。p - 1) + (m + 1)t w ( p - 1) 。1.2 求對(duì)稱矩陣特征值的雅可比法1.2.1 雅可比法求對(duì)稱矩陣特征值的串行算法若矩陣A=aij 是n階實(shí)對(duì)稱矩陣，則存在一

9、個(gè)正交矩陣U，使得UTAU=D其中 D 是一個(gè)對(duì)角矩陣，即él1êD = ê0êMêêë00 Ll2 LMO0L0 ùú0 úM úúln úû這里i (i=1,2,n)是A的特征值，U的第i列是與i 對(duì)應(yīng)的特征向量。雅可比算法主要是通過(guò)正交相似變換將一個(gè)實(shí)對(duì)稱矩陣對(duì)角化，從而求出該矩陣的全部特征值和對(duì)應(yīng)的特征向量。 因此可以用一系列的初等正交變換逐步消去A的非對(duì)角線元素，從而使矩陣A對(duì)角化。設(shè)初 等正交矩陣為R(p,q,)，其(p,p)( q,q)位置的

10、數(shù)據(jù)是cos，(p, q)( q, p)位置的數(shù)據(jù)分別是-sin和 sin(p q)，其它位置的數(shù)據(jù)和一個(gè)同階數(shù)的單位矩陣相同。顯然可以得到：R(p,q,) TR(p,q,)=In不妨記B= R(p,q,)TAR(p,q,)，如果將右端展開(kāi)，則可知矩陣B的元素與矩陣A的元素之 間有如下關(guān)系：bpp = app cos2+aqq sin2+apq sin2bqq = app sin2+aqq cos2-apq sin2bpq = (aqq -app )sin2)/2+apq cos2bqp = bpqbpj = apj cos+aqj sinbqj = -apj sin +aqj cos bip

11、 = aip cos+aiq sinbiq = -aip sin +aiq cos bij = aij其中 1 i, j n且i,j p,q。因?yàn)锳為對(duì)稱矩陣，R為正交矩陣，所以B亦為對(duì)稱矩陣。若要求矩陣B的元素bpq =0，則只需令 (aqq -app )sin2)/2+apq cos2=0，即：tg 2q =-a pq(a qq - a pp ) 2在實(shí)際應(yīng)用時(shí)，考慮到并不需要解出 ，而只需要求出 sin2，sin 和 cos 就可以了。如果限制 值在-/2 < 2 /2，則可令qqm = - a pq ,n = 1 ( a- a pp ) ,w = sgn( n )m容易推出：si

12、n 2q = w ,2sin q =2(1 +w,1 - w 2 )cosq =m 21 - sin 2 q+ n 2利用sin2，sin和cos的值，即得矩陣B的各元素。矩陣A經(jīng)過(guò)旋轉(zhuǎn)變換，選定的非主pq對(duì)角元素apq 及aqp（一般是絕對(duì)值最大的）就被消去，且其主對(duì)角元素的平方和增加了 2a 2 ，而非主對(duì)角元素的平方和減少了 2a2 pq ，矩陣元素總的平方和不變。通過(guò)反復(fù)選取主元素apq ， 并作旋轉(zhuǎn)變換，就逐步將矩陣A變?yōu)閷?duì)角矩陣。在對(duì)稱矩陣中共有(n2-n)/2 個(gè)非主對(duì)角元素要 被消去, 而每消去一個(gè)非主對(duì)角元素需要對(duì) 2n個(gè)元素進(jìn)行旋轉(zhuǎn)變換, 對(duì)一個(gè)元素進(jìn)行旋轉(zhuǎn) 變換需要 2

13、次乘法和 1 次加法，若各取一次乘法運(yùn)算時(shí)間或一次加法運(yùn)算時(shí)間為一個(gè)單位時(shí) 間，則消去一個(gè)非主對(duì)角元素需要 3 個(gè)單位運(yùn)算時(shí)間，所以下述算法 21.3 的一輪計(jì)算時(shí)間 復(fù)雜度為 (n2-n)/2*2n*3=3n2(n-1)=O(n3)。算法 21.3單處理器上雅可比法求對(duì)稱矩陣特征值的算法輸入：矩陣An×n ，輸出：矩陣特征值 EigenvalueBegin(1)while (max >) do (1.1) max=a1,2 (1.2)for i=1 to n dofor j= i+1 to n doif (ai,j) >max) then max =ai,j ，p=i，

14、q=jend if end forend for(1.3)Compute: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=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)for j=1 to n doif (j p) and ( j q) then (i)bp,j= a

15、p,j*co1+ aq,j*si1 (ii)bq,j= -ap,j*si1 + aq,j*co1end if end for(1.5)for i=1 to n doEndif(i p) and (i q) then(i)bi, p= ai,p*co1+ ai,q*si1 (ii)bi, q= - ai,p*si1+ ai,q*co1end if end for(1.6)for i=1 to n do for j=1 to n doai,j=bi,jend for end forend while(2)for i=1 to n doEigenvaluei= ai, iend for1.2.2 雅

16、可比法求對(duì)稱矩陣特征值的并行算法串行雅可比算法逐次尋找非主對(duì)角元絕對(duì)值最大的元素的方法并不適合于并行計(jì)算。因 此，在并行處理中，我們每次并不尋找絕對(duì)值最大的非主對(duì)角元消去，而是按一定的順序?qū)?A 中的所有上三角元素全部消去一遍，這樣的過(guò)程稱為一輪。由于對(duì)稱性，在一輪中，A 的 下三角元素也同時(shí)被消去一遍。經(jīng)過(guò)若干輪，可使 A 的非主對(duì)角元的絕對(duì)值減少，收斂于 一個(gè)對(duì)角方陣。具體算法如下：設(shè)處理器個(gè)數(shù)為p，對(duì)矩陣A按行劃分為 2p塊，每塊含有連續(xù)的m/2 行向量，記 m = én / pù ， 這些行塊依次記為A0 ,A1 , ,A2p-1 ，并將A2i 與A2i+1 存放與

17、標(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.1 p=4 時(shí)的雅可比算法求對(duì)稱矩陣特征值的數(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è)行塊被送至同一處理器時(shí)，令編號(hào)為i的行塊 中的每行元素和編號(hào)為j的行塊中的每行元素配對(duì)，以消去相應(yīng)的非主對(duì)角元素，這樣在

18、所 有的行塊都兩兩配對(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)稱前數(shù)據(jù)塊，后一個(gè)行塊簡(jiǎn)稱后數(shù)據(jù)塊)始終保持不變，將后 數(shù)據(jù)塊發(fā)送給 1 號(hào)處理器，作為 1 號(hào)處理器的前數(shù)據(jù)塊。同時(shí)接收 1 號(hào)處理器發(fā)送的后數(shù)據(jù) 塊作為自己的后數(shù)據(jù)塊。p-1 號(hào)處理器首先將后數(shù)據(jù)塊發(fā)送給編號(hào)為

19、 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) 處理器發(fā)送的數(shù)據(jù)，最后再將緩沖區(qū)中的數(shù)據(jù)發(fā)送給偶數(shù)號(hào)處

20、理器。數(shù)據(jù)塊傳送的具體過(guò)程 描述如下：算法 21.4雅可比法求對(duì)稱矩陣特征值的并行過(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ù)塊end if(2)if (my-rank=p-

21、1) and ( my-rank mod 2 0) then/*處理器 p-1 且其為奇數(shù)*/ (2.1)for i=m/2 to m-1 do/* 將后數(shù)據(jù)塊保存在緩沖區(qū) buffer 中*/for j=0 to n-1 dobufferi-m/2,j=ai,jend for end for(2.2)for i=m/2 to m-1 dobuf i-m/2 =biend for(2.3)for i=0 to m/2-1 do/*將前數(shù)據(jù)塊移至后數(shù)據(jù)塊的位置上*/for j=0 to n-1 doai+m/2,j=ai,jend for end for(2.4)for i=0 to m/2-1

22、 dobi+m/2 =biend for(2.5)接收 p-2 號(hào)處理器發(fā)送的數(shù)據(jù)塊作為自己的前數(shù)據(jù)塊(2.6)將 buffer 中的后數(shù)據(jù)塊發(fā)送給編號(hào)為 p-2 的處理器end if(3)if (my-rank=p-1) and ( my-rank mod 2=0) then/*處理器 p-1 且其為偶數(shù)*/ (3.1)將后數(shù)據(jù)塊發(fā)送給編號(hào)為 p-2 的處理器(3.2)for i=0 to m/2-1 do/*將前數(shù)據(jù)塊移至后數(shù)據(jù)塊的位置上*/for j=0 to n-1 doai+m/2,j=ai,jend for end for(3.3)for i=0 to m/2-1 dobi+m/2

23、 =biend for(3.4)接收 p-2 號(hào)處理器發(fā)送的數(shù)據(jù)塊作為自己的前數(shù)據(jù)塊end if(4)if (my-rank p-1) and ( my-rank 0) then/*其它的處理器*/ (4.1)if (my-rank mod 2=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)for

24、i=0 to m-1 do/* 將前后數(shù)據(jù)塊保存在緩沖區(qū) buffer 中*/for j=0 to n-1 dobufferi,j=ai,jend for end for(vi)for i=0 to m-1 dobufi =biend for(vii)接收編號(hào)為 my_rank-1 的處理器發(fā)送的數(shù)據(jù)塊作為自己的前 數(shù)據(jù)塊Endend if(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 的處理器end

25、 if各處理器并行地對(duì)其局部存儲(chǔ)器中的非主對(duì)角元素aij 進(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,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ì)角元素的最大元

26、max，并廣播給所有處理器以決定 是否進(jìn)行下一輪迭代。具體算法框架描述如下：算法 21.5雅可比法求對(duì)稱矩陣特征值的并行算法輸入：矩陣An×n ，輸出：矩陣 A 的主對(duì)角元素即為 A 的特征值Begin對(duì)所有處理器 my_rank(my_rank=0, p-1)同時(shí)執(zhí)行如下的算法:(a)for i=0 to m-1 dobi =myrank*m+i/* b 記錄處理器中的行塊數(shù)據(jù)在原矩陣 A 中的實(shí)際行號(hào)*/end for(b)while (max>) do/* max 為 A 中所有非對(duì)角元最大的絕對(duì)值*/ (1)for i=my_rank*m to (my_rank+1)*

27、m-2 do/*對(duì)本處理器內(nèi)部所有行兩兩配對(duì)進(jìn)行旋轉(zhuǎn)變換*/for j=i+1 to (my_rank+1)*m-1 do(1.1)r=i mod m ,t=j mod m/*i, j 為進(jìn)行旋轉(zhuǎn)變換行(稱為主行)的 實(shí)際行號(hào), r, t 為它們?cè)趬K內(nèi)的相對(duì)行號(hào)*/(1.2)if (ar,j 0) then/*對(duì)四個(gè)主元素的旋轉(zhuǎn)變換*/ (i)Compute:f=-ar,j, g=( at,j- ar,i)/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)，bp

28、p= 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)for v=0 to n-1 do/*對(duì)兩個(gè)主行其余元素的旋轉(zhuǎn)變換*/if (v i) and ( v j) thenbrv = ar,v*cos1 + at,v*sin1at,v = -ar,v* sin1 + at,v* cos1end ifend for(iii)for v=0 to n-1 doif (v i) and ( v j) thenar,v=brvend if

29、end for(iv)for v=0 to m-1 do/*對(duì)兩個(gè)主列在本處理器內(nèi)的其余元素的旋轉(zhuǎn)變換*/if ( v r) and ( v t) thenbiv = av,i*cos1 + av,j*sin1av,j= - av,i* sin1 + av,j* cos1end if end for(v)for v=0 to m-1doif (v r) and ( v t) thenav,i= bivend if end for(vi)Compute:ar,i=bpp , at,j=bqq , ar,j=bpq , at,i=bqp，/*用 temp1 保存本處理器主行的行號(hào)和旋轉(zhuǎn)參數(shù)*/te

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

31、=temp2(v-1)*4+1，i1=(int)temp2(v-1)*4+2, j1=(int)temp2(v-1)*4+3 (ii )if (s1、c1、 i1、 j1 中有一不為 0) thenif (my-rank current) then for z=0 to m-1 doziz=az,i1*c1 + az,j1*s1a z,j1=- az,i1*s1 + az,j1*c1end forfor z=0 to m-1 doaz,i1= zizend for end ifend if(iii)current=current+1end for end forend for(2)for co

32、unter=1 to 2p-2 do/*進(jìn)行 2p-2 次處理器間的數(shù)據(jù)交換, 并對(duì)交換后處理器中所有行兩兩配對(duì) 進(jìn)行旋轉(zhuǎn)變換*/(2.1)Data_exchange( )/*處理器間的數(shù)據(jù)交換*/(2.2)for i=0 to m/2-1 dofor j=m/2 to m-1 do(i) if (ai,bj 0) then/*對(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)，

33、bpp= ai,bi*cos1*cos1+ aj,bj*sin1*sin1+ai,bj*sin2， bqq= ai,bi* sin1*sin1+aj,bj* cos1*cos1-ai,bj*sin2， bpq=0, bqp=0for v=0 to n-1 do /*對(duì)兩個(gè)主行其余元素的旋轉(zhuǎn)變換*/if (v bi) and ( v bj) thenbrv = ai,v*cos1 + aj,v*sin1aj,v = -ai,v* sin1 + aj,v* cos1end if end forfor v=0 to n-1 doif (v bi) and ( v bj) thenai,v=brven

34、d if end forfor v=0 to m-1 do/*對(duì)本處理器內(nèi)兩個(gè)主列的其余元素旋轉(zhuǎn)變換*/if (v i) and ( v j) thenbiv = av, bi*cos1 + av, bj*sin1av, bj = - av, bi* sin1 + av, bj* cos1end if end forfor v=0 to m-1 doif (v i) and ( v j) thenav, bi=bivend if end forCompute:ai, bi=bpp ,aj, bj=bqq ,ai, bj=bpq ,aj, bi=bqp/*用 temp1 保存本處理器主行的行號(hào)和

35、旋轉(zhuǎn)參數(shù)*/temp10=sin1 ,temp11=cos1,temp12=(float)bi , temp13= (float)bjend forelseCompute: temp10=0,temp11= 0, temp12= 0,temp13= 0end if(ii)將所有處理器 temp1 中的旋轉(zhuǎn)參數(shù)及主行的行號(hào)按處理 器編號(hào)連接起來(lái)并廣播給所有處理器,存于 temp2 中(iii)current=0 (iv)for v=1 to p do/*根據(jù) temp2 中的其它處理器的旋轉(zhuǎn)參數(shù)及主行的行號(hào)對(duì)相關(guān)的列在本處理器的部分進(jìn)行旋轉(zhuǎn)變換*/Compute:s1=temp2(v-1)*4+

36、0,c1=temp2(v-1)*4+1,i1=(int)temp2(v-1)*4+2, j1=(int)temp2(v-1)*4+3if (s1、c1、 i1、 j1 中有一不為 0) then if (my-rank current) thenfor z=0 to m-1 doziz=az,i1*c1 + az,j1*s1a z,j1=- az,i1s1 + az,j1*c1end forfor z=0 to m-1 doaz,i1= zizend for end ifend ifcurrent=current+1end for end forend for(3)Data-exchange(

37、 )/*進(jìn)行一輪中的最后一次處理器間的數(shù)據(jù)交換,使數(shù)據(jù)回到原來(lái)的位置*/ (4)localmax=max/*localmax 為本處理器中非對(duì)角元最大的絕對(duì)值*/(5)for i=0 to m-1 do for j=0 to n-1 doif (m*my-rank+i) j) thenEndif (ai,j> localmax) then localmax=ai,j end if end ifend for end for(6)通過(guò) Allreduce 操作求出所有處理器中 localmax 的最大值為新的 max 值end while在上述算法中, 每個(gè)處理器在一輪迭代中要處理 2 p

38、個(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ì) 2n 個(gè)元素進(jìn)行旋 轉(zhuǎn)變換。而對(duì)一個(gè)元素進(jìn)行旋

39、轉(zhuǎn)變換需要 2 次乘法和 1 次加法，若取一次乘法運(yùn)算時(shí)間或一 次加法運(yùn)算時(shí)間為一個(gè)單位時(shí)間，則其需要 3 個(gè)單位運(yùn)算時(shí)間，即一輪迭代的計(jì)算時(shí)間為T13×2 p ×2n ×m2/43n3/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 = 4t s + m(n + 1)t w (4 p - 2) + n(n - 1) / p + 2t s (p - 1) + 2n(n

40、- 1) / p + 1t w ( p - 1)因此雅可比算法求對(duì)矩陣特征值的一輪并行計(jì)算時(shí)間為Tp = T1 + T2 。我們的大量實(shí)驗(yàn)結(jié)果說(shuō)明，對(duì)于n階方陣，用上述算法進(jìn)行并行計(jì)算，一般需要不超過(guò)O(logn)輪就可以收斂。MPI 源程序請(qǐng)參見(jiàn)章末附錄。1.3 求對(duì)稱矩陣特征值的單側(cè)旋轉(zhuǎn)法1.3.1 單側(cè)旋轉(zhuǎn)法的算法描述求解對(duì)稱方陣特征值及特征向量的雅可比算法在每次消去計(jì)算前都要對(duì)非主對(duì)角元素 選擇最大元，導(dǎo)致非實(shí)際計(jì)算開(kāi)銷很大。在消去計(jì)算時(shí)，必須對(duì)方陣同時(shí)進(jìn)行行、列旋轉(zhuǎn)變 換，這稱之為雙側(cè)旋轉(zhuǎn)(Two-side Rotation)變換。在雙側(cè)旋轉(zhuǎn)變換中，方陣的行、列方向都 有數(shù)據(jù)相關(guān)關(guān)系

41、，使得整個(gè)計(jì)算中的相關(guān)關(guān)系復(fù)雜，特別是在對(duì)高階矩陣進(jìn)行特征值求解時(shí)， 增加了處理器間的通信開(kāi)銷。針對(duì)傳統(tǒng)雅可比算法的缺點(diǎn)，可以將雙側(cè)旋轉(zhuǎn)改為單側(cè)旋轉(zhuǎn)， 得出一種求對(duì)稱矩陣特征值的快速算法。這一算法對(duì)矩陣僅實(shí)施列變換，使得數(shù)據(jù)相關(guān)關(guān)系 僅在同列之間，因此方便數(shù)據(jù)劃分，適合并行計(jì)算，稱為單側(cè)旋轉(zhuǎn)法（One-side Rotation）。 若A 為一對(duì)稱方陣， 是A 的特征值，x 是A 的特征向量，則有：Ax=x 。又A=AT ，所以 ATAx=Ax=x，這說(shuō)明了 2是ATA的特征值，x是ATA的特征向量 ，即ATA的特征值是A的特 征值的平方，并且它們的特征向量相同。我們使用 18.7.1 節(jié)中所

42、介紹的Givens旋轉(zhuǎn)變換對(duì)A進(jìn)行一系列的列變換，得到方陣Q使 其各列兩兩正交，即AVQ，這里V為表示正交化過(guò)程的變換方陣。因QTQ=(AV) TAV=V TATAV= diag(1 ,2 , ,n ) 為n階對(duì)角方陣，可見(jiàn)這里1 ,2 , ,n 是矩陣ATA的特征值，V的各列是ATA的特征向量。由此可得：d = l 2(i=1,2, ,n)。設(shè)Q的第i列為q , 則 q Tq = ,iiiiiii|qi |2 =qT qi =d i = li。因此在將A進(jìn)行列變換得到各列兩兩正交的方陣Q后，其各列的譜范數(shù)|qi |2 即為A的特征值的絕對(duì)值。記V的第i列為vi ， AVQ，所以Avi = q

43、i 。又因?yàn)関i 為 A的特征向量，所以Avi =i vi ，即推得i vi = qi 。因此i 的符號(hào)可由向量 qi 及vi 的對(duì)應(yīng)分量是 否同號(hào)判別，實(shí)際上在具體算法中我們只要判斷它們的第一個(gè)分量是否同號(hào)即可。若相同， 則i 取正值，否則取負(fù)值。求對(duì)稱矩陣特征值的單側(cè)旋轉(zhuǎn)法的串行算法如下： 算法 21.6求對(duì)稱矩陣特征值的單側(cè)旋轉(zhuǎn)法輸入：對(duì)稱矩陣 A，精度 e輸出：矩陣 A 的特征值存于向量 B 中Begin(1)while (p>)p=0for i=1 to n dofor j=i+1 to n do(1.1 )sum0=0,sum1=0,sum2=0 (1.2 )for k=1

44、to n dosum0= sum0+ak,i*ak, j sum1= sum1+ ak,i* ak,i sum2= sum2+ ak, j* ak, jend for(1.3 )if (sum0>) then(i) aa=2*sum0 bb=sum1-sum2 rr=sqrt(aa*aa+bb*bb)(ii) if (bb0) thenc=sqrt(bb+rr)/(2*rr)s=aa/(2*rr*c)elseend ifs=sqrt(rr-bb)/(2*rr)c=aa/(2*rr*s)(iii) for k=1 to n do tempk=c*ak,i+s*ak,j ak,j=-s*ak

45、,i+c*ak,jend for(iv) for k=1 to n dotemp1k= c*ek,i+s*ek,jek,j=-s*ek,i+c*ek,jend for(v) for k=1 to n doend for end forend while(2)for i=1 to n dosu=0ak,i= tempkend for(vi) for k=1 to n doek,i= temp1kend for(vii) if (sum0>p) then p=sum0 end if end ifEndfor j=1 to n dosu=su+aj,i* aj,iend forBj=sqrts

46、uif (e0,j*a0,j<0) then Bj= - Bjendif end for上述算法的一輪迭代需進(jìn)行n(n-1)/2 次旋轉(zhuǎn)變換，若取一次乘法或一次加法運(yùn)算時(shí)間為 一個(gè)單位時(shí)間，則一次旋轉(zhuǎn)變換要做 3 次內(nèi)積運(yùn)算而消耗 6n個(gè)單位時(shí)間；與此同時(shí)，兩列 元素進(jìn)行正交計(jì)算還需要 12n 個(gè)單位時(shí)間，所以奇異值分解的一輪計(jì)算時(shí)間復(fù)雜度為 n(n-1)/2*(12 n +6n)= 9n(n-1) n=O(n3)。1.3.2 求對(duì)稱矩陣特征值的單側(cè)旋轉(zhuǎn)法的并行計(jì)算在求對(duì)稱矩陣特征值的單側(cè)旋轉(zhuǎn)法的計(jì)算中, 主要的計(jì)算是矩陣的各列正交化過(guò)程.為 了進(jìn)行并行計(jì)算, 我們對(duì) n 階對(duì)稱矩陣 A

47、 按行劃分為 p 塊(p 為處理器數(shù))，每塊含有連續(xù)的 q 行向量，這里 q=n/p，第 i 塊包含 A 的第 i×q, ,(i+1)×q-1 行向量，其數(shù)據(jù)元素被分配到 第 i 號(hào)處理器上(i=0,1,p-1)。E 矩陣取階數(shù)為 n 的單位矩陣 I，按同樣方式進(jìn)行數(shù)據(jù)劃分， 每塊含有連續(xù)的 q 行向量。對(duì)矩陣 A 的每一個(gè)列向量而言，它被分成 p 個(gè)長(zhǎng)度為 q 的子向 量，分布于 p 個(gè)處理器之中，它們協(xié)同地對(duì)各列向量做正交計(jì)算。在對(duì)第 i 列與第 j 列進(jìn)行 正交計(jì)算時(shí)，各個(gè)處理器首先求其局部存儲(chǔ)器中的 q 維子向量i,j、i,i、j,j的內(nèi)積，然后 通過(guò)歸約操作的求和

48、運(yùn)算求得整個(gè) n 維向量對(duì)i,j、i,i、j,j的內(nèi)積并廣播給所有處理器， 最后各處理器利用這些內(nèi)積對(duì)其局部存儲(chǔ)器中的第 i 列及第 j 列 q 維子向量的元素做正交計(jì) 算。算法 21.7 按這樣的方式對(duì)所有列對(duì)正交化一次以完成一輪運(yùn)算，重復(fù)進(jìn)行若干輪運(yùn)算， 直到迭代收斂為止。在各列正交后, 編號(hào)為 0 的處理器收集各處理器中的計(jì)算結(jié)果，由 0 號(hào) 處理器計(jì)算矩陣的特征值. 具體算法框架描述如下：算法 21.7求對(duì)稱矩陣特征值的并行單側(cè)旋轉(zhuǎn)算法 輸入：對(duì)稱矩陣 A，精度 e輸出：對(duì)稱矩陣 A 的特征值存于向量 B 中Begin對(duì)所有處理器 my_rank(my_rank=0, p-1)同時(shí)執(zhí)行

49、如下的算法：(a) while (not convergence) do(1)k=0(2)for i=0 to n-1 do(2.1)for j=i+1 to n-1 do(i)sum0=0,sum1=0,sum2=0 (ii)計(jì)算本處理器所存的列子向量 a*, i、a*, j的內(nèi)積賦值給 ss0 (iii)計(jì)算本處理器所存的列子向量 a*, i、a*, i的內(nèi)積賦值給 ss1 (iv)計(jì)算本處理器所存的列子向量 a*, j、a*, j的內(nèi)積賦值給 ss2 (v)通過(guò)規(guī)約操作實(shí)現(xiàn):計(jì)算所有處理器中 ss0的和賦給 sum 0計(jì)算所有處理器中 ss1的和賦給 sum 1計(jì)算所有處理器中 ss2的

50、和賦給 sum2將 sum 0、sum 1、sum 2 的值廣播到所有處理器中(vi)if (sum(0)>e ) then/*各個(gè)處理器并行進(jìn)行對(duì) i 和 j 列正交化*/aa=2*sum0bb=sum1-sum2rr=sqrt(aa*aa+bb*bb)if (bb>=0) then c=sqrt(bb+rr)/(2*rr) s=aa/(2*rr*c)elses=sqrt(rr-bb)/(2*rr)c=aa/(2*rr*s)end iffor v=0 to q-1 do tempv=c*av, i+s*av, j av, j=(-s)*av, i+c*av, jend forfo

51、r v=0 to z-1 dotemp1v=c*ev, i+s*ev, j ev,j=(-s)*ev, i+c*ev, jend forfor v=0 to q-1 doav,i=tempvend forfor v=0 to z-1 doev,i=temp1vend fork k+1end if end forend for end while0 號(hào)處理器從其余各處理器中接收子矩陣 a 和 e，得到經(jīng)過(guò)變換的矩陣 A 和 I：/*0 號(hào)處理器負(fù)責(zé)計(jì)算特征值*/ (b) for i=1 to n doEnd(1)su=0(2)for j=1 to n dosu=su+Aj,i* Aj,iend

52、for(3)Bj=sqrt(su)(4)if (I0,j*a0,j<0) then Bj= - Bjendif end for上述算法的一輪迭代要進(jìn)行n(n-1)/2 次旋轉(zhuǎn)變換, 而一次旋轉(zhuǎn)變換需要做 3 次內(nèi)積運(yùn)算， 若取一次乘法或一次加法運(yùn)算時(shí)間為一個(gè)單位時(shí)間，則需要 6q個(gè)單位時(shí)間，另外，還要對(duì) 四列子向量中元素進(jìn)行正交計(jì)算花費(fèi) 12q 個(gè)單位時(shí)間，所以一輪迭代需要的計(jì)算時(shí)間 T1 =n(n-1)6q；內(nèi)積計(jì)算需要通信，一輪迭代共做歸約操作n(n-1)/2 次，每次通信量為 3，因而通信時(shí)間T2 = 2t s (Tp =T1 +T2 。p - 1) + 3t w ( p - 1) * n (n - 1)2 。由此得出一輪迭代的并行計(jì)算時(shí)間MPI 源程序請(qǐng)參見(jiàn)所附光盤。1.4 求一般矩陣全部特征值的QR方法1.4.1 QR方法求一般矩陣全部特征值的串行算法在 18.6 節(jié)中，我們介紹了對(duì)一個(gè)n階方陣A