并行計算矩陣特征值計算

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

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

3、若vk-1 與vk 的各個分量同號且成比例，則1 =uk ；若vk-1 與vk 的各個分量異號且成比 例，則1 = -uk 。若各取一次乘法和加法運算時間、一次除法運算時間、一次比較運算時 間為一個單位時間，則因為一輪計算要做一次矩陣向量相乘、一次求最大元操作和一次規(guī)格 化操作，所以下述乘冪法串行算法 21.1 的一輪計算時間為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 乘冪法并行算法乘冪法求矩陣特征值由反復進行矩陣向量相乘來實現(xiàn)，因而可以采用矩陣向量相乘的數(shù) 據(jù)劃分方法。設處理器個數(shù)為 p，對矩陣 A 按行劃分為 p 塊，

5、每塊含有連續(xù)的 m 行向量， 這里 m = én / pù ，編號為 i 的處理器含有 A 的第 im 至第(i+1)m-1 行數(shù)據(jù)，(i=0,1, ,p-1)，初 始向量 v 被廣播給所有處理器。各處理器并行地對存于局部存儲器中 a 的行塊和向量 v 做乘積操作，并求其 m 維結果 向量中的最大值 localmax，然后通過歸約操作的求最大值運算得到整個 n 維結果向量中的最 大值 max 并廣播給所有處理器，各處理器利用 max 進行規(guī)格化操作。最后通過擴展收集操 作將各處理器中的 m 維結果向量按處理器編號連接起來并廣播給所有處理器，以進行下一 次迭代計算，直至收斂。

6、具體算法框架描述如下：算法 21.2乘冪法求解矩陣最大特征值的并行算法輸入：系數(shù)矩陣An×n ，初始向量v n×1 ，輸出：最大的特征值 maxBegin對所有處理器 my_rank(my_rank=0, p-1)同時執(zhí)行如下的算法:while (diff>) do/* diff 為特征向量的各個分量新舊值之差的最大值*/ (1)for i=0 to m-1 do/*對所存的行計算特征向量的相應分量*/(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/*對所計算的特征向量的相應分量 求新舊值之差的最大值 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/*對所計算的特征向量歸一化 */bi =bi/maxend for(6)用 Allgather 操作將各個處理器中計算出的特征向量的分量的新值組合并廣播到 所有處理器中(7)diff=max-oldmax, oldmax=m

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

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

10、數(shù)據(jù)是cos，(p, q)( q, p)位置的數(shù)據(jù)分別是-sin和 sin(p q)，其它位置的數(shù)據(jù)和一個同階數(shù)的單位矩陣相同。顯然可以得到：R(p,q,) TR(p,q,)=In不妨記B= R(p,q,)TAR(p,q,)，如果將右端展開，則可知矩陣B的元素與矩陣A的元素之 間有如下關系：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。因為A為對稱矩陣，R為正交矩陣，所以B亦為對稱矩陣。若要求矩陣B的元素bpq =0，則只需令 (aqq -app )sin2)/2+apq cos2=0，即：tg 2q =-a pq(a qq - a pp ) 2在實際應用時，考慮到并不需要解出 ，而只需要求出 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)過旋轉變換，選定的非主pq對角元素apq 及aqp（一般是絕對值最大的）就被消去，且其主對角元素的平方和增加了 2a 2 ，而非主對角元素的平方和減少了 2a2 pq ，矩陣元素總的平方和不變。通過反復選取主元素apq ， 并作旋轉變換，就逐步將矩陣A變?yōu)閷蔷仃嚒Ｔ趯ΨQ矩陣中共有(n2-n)/2 個非主對角元素要 被消去, 而每消去一個非主對角元素需要對 2n個元素進行旋轉變換, 對一個元素進行旋轉 變換需要 2

13、次乘法和 1 次加法，若各取一次乘法運算時間或一次加法運算時間為一個單位時 間，則消去一個非主對角元素需要 3 個單位運算時間，所以下述算法 21.3 的一輪計算時間 復雜度為 (n2-n)/2*2n*3=3n2(n-1)=O(n3)。算法 21.3單處理器上雅可比法求對稱矩陣特征值的算法輸入：矩陣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、可比法求對稱矩陣特征值的并行算法串行雅可比算法逐次尋找非主對角元絕對值最大的元素的方法并不適合于并行計算。因 此，在并行處理中，我們每次并不尋找絕對值最大的非主對角元消去，而是按一定的順序將 A 中的所有上三角元素全部消去一遍，這樣的過程稱為一輪。由于對稱性，在一輪中，A 的 下三角元素也同時被消去一遍。經(jīng)過若干輪，可使 A 的非主對角元的絕對值減少，收斂于 一個對角方陣。具體算法如下：設處理器個數(shù)為p，對矩陣A按行劃分為 2p塊，每塊含有連續(xù)的m/2 行向量，記 m = én / pù ， 這些行塊依次記為A0 ,A1 , ,A2p-1 ，并將A2i 與A2i+1 存放與

17、標號為i的處理器中。每輪計算開始，各處理器首先對其局部存儲器中所存的兩個行塊的所有行兩兩配對進行 旋轉變換，消去相應的非對角線元素。然后按圖 21.1 所示規(guī)律將數(shù)據(jù)塊在不同處理器之間傳送，以消去其它非主對角元素。圖 1.1 p=4 時的雅可比算法求對稱矩陣特征值的數(shù)據(jù)交換圖這里長方形框中兩個方格內(nèi)的整數(shù)被看成是所移動的行塊的編號。在要構成新的行塊配 對時，只要將每個處理器所對應的行塊按箭頭方向移至相鄰的處理器即可，這樣的傳送可以 在行塊之間實現(xiàn)完全配對。當編號為i和j的兩個行塊被送至同一處理器時，令編號為i的行塊 中的每行元素和編號為j的行塊中的每行元素配對，以消去相應的非主對角元素，這樣在

18、所 有的行塊都兩兩配對之后，可以將所有的非主對角元素消去一遍，從而完成一輪計算。由圖1.1可以看出，在一輪計算中，處理器之間要 2p-1 次交換行塊。為保證不同行塊配對時可以 將原矩陣A的非主對角元素消去，引入變量b記錄每個處理器中的行塊數(shù)據(jù)在原矩陣A中的實 際行號。在行塊交換時，變量b也跟著相應的行塊在各處理器之間傳送。處理器之間的數(shù)據(jù)塊交換存在如下規(guī)律：0 號處理器前一個行塊(簡稱前數(shù)據(jù)塊，后一個行塊簡稱后數(shù)據(jù)塊)始終保持不變，將后 數(shù)據(jù)塊發(fā)送給 1 號處理器，作為 1 號處理器的前數(shù)據(jù)塊。同時接收 1 號處理器發(fā)送的后數(shù)據(jù) 塊作為自己的后數(shù)據(jù)塊。p-1 號處理器首先將后數(shù)據(jù)塊發(fā)送給編號為

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

20、理器。數(shù)據(jù)塊傳送的具體過程 描述如下：算法 21.4雅可比法求對稱矩陣特征值的并行過程中處理器之間的數(shù)據(jù)塊交換算法 輸入：矩陣 A 的行數(shù)據(jù)塊和向量 b 的數(shù)據(jù)塊分布于各個處理器中 輸出：在處理器陣列中傳送后的矩陣 A 的行數(shù)據(jù)塊和向量 b 的數(shù)據(jù)塊Begin對所有處理器 my_rank(my_rank=0, p-1)同時執(zhí)行如下的算法:/*矩陣 A 和向量 b 為要傳送的數(shù)據(jù)塊*/ (1)if (my-rank=0) then/*0 號處理器*/(1.1)將后數(shù)據(jù)塊發(fā)送給 1 號處理器(1.2)接收 1 號處理器發(fā)送來的后數(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 號處理器發(fā)送的數(shù)據(jù)塊作為自己的前數(shù)據(jù)塊(2.6)將 buffer 中的后數(shù)據(jù)塊發(fā)送給編號為 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ā)送給編號為 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 號處理器發(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ù)號處理器*/(i)將前數(shù)據(jù)塊發(fā)送給編號為 my_rank+1 的處理器(ii)將后數(shù)據(jù)塊發(fā)送給編號為 my_rank-1 的處理器(ii)接收編號為 my_rank-1 的處理器發(fā)送的數(shù)據(jù)塊作為自己的前 數(shù)據(jù)塊(iv)接收編號為 my_rank+1 的處理器發(fā)送的數(shù)據(jù)塊作為自己的后 數(shù)據(jù)塊else/*奇數(shù)號處理器*/(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)接收編號為 my_rank-1 的處理器發(fā)送的數(shù)據(jù)塊作為自己的前 數(shù)據(jù)塊Endend if(viii)接收編號為 my_rank+1 的處理器發(fā)送的數(shù)據(jù)塊作為自己的 后數(shù)據(jù)塊(ix)將存于 buffer 中的前數(shù)據(jù)塊發(fā)送給編號為 my_rank+1 的處理 器(x)將存于 buffer 中的后數(shù)據(jù)塊發(fā)送給編號為 my_rank-1 的處理器end

25、 if各處理器并行地對其局部存儲器中的非主對角元素aij 進行消去，首先計算旋轉參數(shù)并對 第i行和第j行兩行元素進行旋轉行變換。然后通過擴展收集操作將相應的旋轉參數(shù)及第i列和 第j列的列號按處理器編號連接起來并廣播給所有處理器。各處理器在收到這些旋轉參數(shù)和 列號之后，按 0,1,p-1 的順序依次讀取旋轉參數(shù)及列號并對其m行中的第i列和第j列元素進 行旋轉列變換。經(jīng)過一輪計算的 2p-1 次數(shù)據(jù)交換之后，原矩陣 A 的所有非主對角元素都被消去一次。此時，各處理器求其局部存儲器中的非主對角元素的最大元 localmax，然后通過歸約操作的 求最大值運算求得將整個 n 階矩陣非主對角元素的最大元

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

27、m-2 do/*對本處理器內(nèi)部所有行兩兩配對進行旋轉變換*/for j=i+1 to (my_rank+1)*m-1 do(1.1)r=i mod m ,t=j mod m/*i, j 為進行旋轉變換行(稱為主行)的 實際行號, r, t 為它們在塊內(nèi)的相對行號*/(1.2)if (ar,j 0) then/*對四個主元素的旋轉變換*/ (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/*對兩個主行其余元素的旋轉變換*/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/*對兩個主列在本處理器內(nèi)的其余元素的旋轉變換*/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 保存本處理器主行的行號和旋轉參數(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 中的旋轉參數(shù)及主行的行號 按處理器編號連接起來并廣播給所有處理器,存于 temp2 中(1.4)current=0(1.5) for v=1 to p do/*根據(jù) temp2 中的其它處理器的旋轉參數(shù)及主行的行號對相關的 列在本處理器的部分進行旋轉變換*/(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/*進行 2p-2 次處理器間的數(shù)據(jù)交換, 并對交換后處理器中所有行兩兩配對 進行旋轉變換*/(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/*對四個主元素的旋轉變換*/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 /*對兩個主行其余元素的旋轉變換*/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/*對本處理器內(nèi)兩個主列的其余元素旋轉變換*/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 保存本處理器主行的行號和

35、旋轉參數(shù)*/temp10=sin1 ,temp11=cos1,temp12=(float)bi , temp13= (float)bjend forelseCompute: temp10=0,temp11= 0, temp12= 0,temp13= 0end if(ii)將所有處理器 temp1 中的旋轉參數(shù)及主行的行號按處理 器編號連接起來并廣播給所有處理器,存于 temp2 中(iii)current=0 (iv)for v=1 to p do/*根據(jù) temp2 中的其它處理器的旋轉參數(shù)及主行的行號對相關的列在本處理器的部分進行旋轉變換*/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、 )/*進行一輪中的最后一次處理器間的數(shù)據(jù)交換,使數(shù)據(jù)回到原來的位置*/ (4)localmax=max/*localmax 為本處理器中非對角元最大的絕對值*/(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)通過 Allreduce 操作求出所有處理器中 localmax 的最大值為新的 max 值end while在上述算法中, 每個處理器在一輪迭代中要處理 2 p

38、個行塊對, 由于每一個行塊對含有 m/2 行, 因而對每一個行塊對的處理要有(m/2)2=m2/4 個行的配對, 即消去m2/4 個非主對角 元素. 每消去一個非主對角元素要對同行n 個元素和同列n 個元素進行旋轉變換. 由于按行 劃分數(shù)據(jù)塊, 對同行n 個元素進行旋轉變換的過程是在本處理器中進行的. 而對同列n 個元 素進行旋轉變換的過程則分布在所有處理器中進行. 但由于所有處理器同時在為其它處理 器的元素消去過程進行列的旋轉變換，對每個處理器而言，對列元素進行旋轉變換的處理總 量仍然是n 個元素。 因此，一個處理器每消去一個非主對角元素共要對 2n 個元素進行旋 轉變換。而對一個元素進行旋

39、轉變換需要 2 次乘法和 1 次加法，若取一次乘法運算時間或一 次加法運算時間為一個單位時間，則其需要 3 個單位運算時間，即一輪迭代的計算時間為T13×2 p ×2n ×m2/43n3/p；在每輪迭代中，各個處理器之間以點對點的通信方式相互錯開交換數(shù)據(jù) 2p-1 次，通信量為mn+m，擴展收集操作n(n-1)/(2p)次，通信量為 4，另外有歸 約操作 1 次，通信量為 1 ，從而不難得出上述算法求解過程中的總通信時間為：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)因此雅可比算法求對矩陣特征值的一輪并行計算時間為Tp = T1 + T2 。我們的大量實驗結果說明，對于n階方陣，用上述算法進行并行計算，一般需要不超過O(logn)輪就可以收斂。MPI 源程序請參見章末附錄。1.3 求對稱矩陣特征值的單側旋轉法1.3.1 單側旋轉法的算法描述求解對稱方陣特征值及特征向量的雅可比算法在每次消去計算前都要對非主對角元素 選擇最大元，導致非實際計算開銷很大。在消去計算時，必須對方陣同時進行行、列旋轉變 換，這稱之為雙側旋轉(Two-side Rotation)變換。在雙側旋轉變換中，方陣的行、列方向都 有數(shù)據(jù)相關關系

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

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

43、i 。又因為vi 為 A的特征向量，所以Avi =i vi ，即推得i vi = qi 。因此i 的符號可由向量 qi 及vi 的對應分量是 否同號判別，實際上在具體算法中我們只要判斷它們的第一個分量是否同號即可。若相同， 則i 取正值，否則取負值。求對稱矩陣特征值的單側旋轉法的串行算法如下： 算法 21.6求對稱矩陣特征值的單側旋轉法輸入：對稱矩陣 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上述算法的一輪迭代需進行n(n-1)/2 次旋轉變換，若取一次乘法或一次加法運算時間為 一個單位時間，則一次旋轉變換要做 3 次內(nèi)積運算而消耗 6n個單位時間；與此同時，兩列 元素進行正交計算還需要 12n 個單位時間，所以奇異值分解的一輪計算時間復雜度為 n(n-1)/2*(12 n +6n)= 9n(n-1) n=O(n3)。1.3.2 求對稱矩陣特征值的單側旋轉法的并行計算在求對稱矩陣特征值的單側旋轉法的計算中, 主要的計算是矩陣的各列正交化過程.為 了進行并行計算, 我們對 n 階對稱矩陣 A

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

48、運算求得整個 n 維向量對i,j、i,i、j,j的內(nèi)積并廣播給所有處理器， 最后各處理器利用這些內(nèi)積對其局部存儲器中的第 i 列及第 j 列 q 維子向量的元素做正交計 算。算法 21.7 按這樣的方式對所有列對正交化一次以完成一輪運算，重復進行若干輪運算， 直到迭代收斂為止。在各列正交后, 編號為 0 的處理器收集各處理器中的計算結果，由 0 號 處理器計算矩陣的特征值. 具體算法框架描述如下：算法 21.7求對稱矩陣特征值的并行單側旋轉算法 輸入：對稱矩陣 A，精度 e輸出：對稱矩陣 A 的特征值存于向量 B 中Begin對所有處理器 my_rank(my_rank=0, p-1)同時執(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)計算本處理器所存的列子向量 a*, i、a*, j的內(nèi)積賦值給 ss0 (iii)計算本處理器所存的列子向量 a*, i、a*, i的內(nèi)積賦值給 ss1 (iv)計算本處理器所存的列子向量 a*, j、a*, j的內(nèi)積賦值給 ss2 (v)通過規(guī)約操作實現(xiàn):計算所有處理器中 ss0的和賦給 sum 0計算所有處理器中 ss1的和賦給 sum 1計算所有處理器中 ss2的

50、和賦給 sum2將 sum 0、sum 1、sum 2 的值廣播到所有處理器中(vi)if (sum(0)>e ) then/*各個處理器并行進行對 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 號處理器從其余各處理器中接收子矩陣 a 和 e，得到經(jīng)過變換的矩陣 A 和 I：/*0 號處理器負責計算特征值*/ (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上述算法的一輪迭代要進行n(n-1)/2 次旋轉變換, 而一次旋轉變換需要做 3 次內(nèi)積運算， 若取一次乘法或一次加法運算時間為一個單位時間，則需要 6q個單位時間，另外，還要對 四列子向量中元素進行正交計算花費 12q 個單位時間，所以一輪迭代需要的計算時間 T1 =n(n-1)6q；內(nèi)積計算需要通信，一輪迭代共做歸約操作n(n-1)/2 次，每次通信量為 3，因而通信時間T2 = 2t s (Tp =T1 +T2 。p - 1) + 3t w ( p - 1) * n (n - 1)2 。由此得出一輪迭代的并行計算時間MPI 源程序請參見所附光盤。1.4 求一般矩陣全部特征值的QR方法1.4.1 QR方法求一般矩陣全部特征值的串行算法在 18.6 節(jié)中，我們介紹了對一個n階方陣A