并行计算-矩阵特征值计算--

1、9矩阵特征值计算在实际的工程计算中，经常会遇到求n阶方阵A的特征值(Eigenvalue)与特征向量(Eigenvector)的问题。对于一个方阵A,如果数值入使方程组Ax=入x即(A-Xn)x=0有非零解向量(SolutionVector)x,则称讷方阵A的特征值，而非零向量x为特征值所对应的特征向量，其中In为n阶单位矩阵。由于根据定义直接求矩阵特征值的过程比较复杂，因此在实际计算中，往往采取一些数值方法。本章主要介绍求一般方阵绝对值最大的特征值的乘哥(Power)法、求对称方阵特征值的雅可比法和单侧旋转(One-sideRotation)法以及求一般矩阵全部特征值的QR方法及一些相关的并

2、行算法。1.1 求解矩阵最大特征值的乘哥法1.1.1 乘幕法及其用行算法在许多实际问题中，只需要计算绝对值最大的特征值，而并不需要求矩阵的全部特征值。乘哥法是一种求矩阵绝对值最大的特征值的方法。记实方阵A的n个特征值为%i=(1,2,n),且满足：1213-n特征值入对应的特征向量为xi。乘哥法的做法是：取n维非零向量V0作为初始向量；对于k=1,2,做如下迭代：直至uk1Uk.,|Uk=AVk-1Vk=Uk/UkJ00卜II的止，这时Vk+i就是A的绝对值最大的特征值入i所对应的特征向8量Xi。若Vk-1与Vk的各个分量同号且成比例，则入产IUkI若Vk-1与Vk的各个分量异号且成比例,则入

3、产-IUkIooo若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时间，则因为一轮计算要做一次矩阵向量相乘、一次求最大元操作和一次规格化操作，所以下述乘哥法串行算法21.1的一轮计算时间为n2+2n=O(n2)。算法21.1单处理器上乘哥法求解矩阵最大特征值的算法输入：系数矩阵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乘幕法并行算法乘哥法求矩阵特征值由反复进行矩阵向量相乘来实现，因而可以采用矩阵向量相乘的数据划分方法。设处理器个数为p,对矩阵A按行划分为p块，每块含有连续的m行向量，这里m=n/p,编号为i的处理器含有A的第im至第(i+1)m-1行数据，(i=0,1,p-1),初始向量v被广播给所有处理器。各处理器并行地对存于局部存储器中a的行块和向量v做乘积操作，并求其m维结果向量中的最大值lo

5、calmax,然后通过归约操作的求最大值运算得到整个n维结果向量中的最大值max并广播给所有处理器，各处理器利用max进行规格化操作。最后通过扩展收集操作将各处理器中的m维结果向量按处理器编号连接起来并广播给所有处理器，以进行下一次迭代计算，直至收敛。具体算法框架描述如下：算法21.2乘哥法求解矩阵最大特征值的并行算法输入：系数矩阵An。，初始向量Vnx1,e输出：最大的特征值maxBegin对所有处理器my_rank(my_rank=0,，p-1)同时执行如下的算法：while(diff)sdo/*diff为特征向量的各个分量新旧值之差的最大值*/(1)fori=0tom-1do/*对所存的

6、行计算特征向量的相应分量*/(1.(1) m=0(1.(2) forj=0ton-1dosum=sum+ai,j*xjendfor(1.(3) i=sumendfor(2)localmax=b0/*对所计算的特征向量的相应分量求新旧值之差的最大值localmax*/(3)fori=1tom-1doif(bi：localmax)thenlocalmax=biendifendfor(4)用Allreduce操作求出所有处理器中localmax值的最大值max并广播到所有处理器中(5)fori=0tom-1do/*对所计算的特征向量归一化*/bi=bi/maxendfor(6)用Allgather操

7、作将各个处理器中计算出的特征向量的分量的新值组合并广播到所有处理器中(7)diff=max-oldmax,oldmax=maxendwhileEnd若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时间，一轮迭代的计算时间为mn+2m;一轮迭代中，各处理器做一次归约操作，通信量为1,一次扩展收集操作，通信量为m,则通信时间为4ts(pp-1)十(m+1)tw(p-1)o因此乘哥法的一轮并行计算时间为Tp=mn+2m+4ts/-p-1)+(m+1)tw(p-1)。(、MPI源程序请参见所附光盘。1.2求对称矩阵特征值的雅可比法1.2.1 雅可比法求对称矩阵特征值的串行算法

8、若矩阵A=aj是n阶实对称矩阵，则存在一个正交矩阵U,使得UTAU=D其中D是一个对角矩阵，即入o0D=0屹0I0?九n0这里入(i=1,2,n)是A的特征值，U的第i列是与入对应的特征向量。雅可比算法主要是通过正交相似变换将一个实对称矩阵对角化，从而求出该矩阵的全部特征值和对应的特征向量。因此可以用一系列的初等正交变换逐步消去A的非对角线元素，从而使矩阵A对角化。设初等正交矩阵为R(p,q,其(p,p)(q,q)位置的数据是cos。，(p,q)(q,p)位置的数据分别是-sin酥口sin&pwq),其它位置的数据和一个同阶数的单位矩阵相同。显然可以得到：R(p,q,)TR(p,q,0)=In

9、不妨记8=R(p,q,9)TAR(p,q,0),如果将右端展开，则可知矩阵B的元素与矩阵A的元素之间有如下关系：bpp=appcos2卅aqqsin20apqsin20bpq=(aqq-app)sin20)/2+apqcos20bpj=apjcos卅aqjsin0bip=aipcos卅aiqsin022八bqq一appsin9+aqqcos9-apqsin20bqp=bpqbqj=-apjsin0+aqjcos0biq=-aipsin0+aiqcos0bij=aijR为正交矩阵，所以B亦为对称矩阵。若其中1wi,jwn且i,jwp,q。因为A为对称矩阵,要求上I阵B的元素bpq=0,则只需令(

10、aqq-app)sin20)/2+apqcos29=0,即:pqtg2u=(aqq-app)2sin20,sin0和cos0就可以了。在实际应用时，考虑到并不需要解出&而只需要求出如果限制。值在-兀/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雅可比法求对称矩阵特征值的并行算法串行雅可比算法逐次寻找非主对角元绝对值最大的元素的方法并不适合于并行计算。此，在并行处理中，我们每次并不寻找绝对值最大的非主对角元消去，而是按一定的顺序将中的所有上三角元素全部消去一遍，这样的过程称为一轮。由于对称性，在一轮中，A的三角元素也同时被消去一遍。经过若干轮，可使A的非主对角元的绝对值减少，收敛于个对角方阵。具体算法如下：设处理器

13、个数为p对矩阵A按行划分为2P块每块含有连续的m/2行向量，记m=n/p,些行块依次记为Ao,Ai,A2P，并将A2i与A2i+1存放与标号为i的处理器中。每轮计算开始，各处理器首先对其局部存储器中所存的两个行块的所有行两两配对进行旋转变换，消去相应的非对角线元素。然后按图21.1所示规律将数据块在不同处理器之间传送，以消去其它非主对角元素。图1.1p=4时的雅可比算法求对称矩阵特征值的数据交换图这里长方形框中两个方格内的整数被看成是所移动的行块的编号。在要构成新的行块配对时，只要将每个处理器所对应的行块按箭头方向移至相邻的处理器即可，这样的传送可以在行块之间实现完全配对。当编号为i和j的两个

14、行块被送至同一处理器时，令编号为i的行块中的每行元素和编号为j的行块中的每行元素配对，以消去相应的非主对角元素，这样在所有的行块都两两配对之后，可以将所有的非主对角元素消去一遍，从而完成一轮计算。由图1.1可以看出，在一轮计算中，处理器之间要2p-1次交换行块。为保证不同行块配对时可以将原矩阵A的非主对角元素消去，引入变量b记录每个处理器中的行块数据在原矩阵A中的实际行号。在行块交换时，变量b也跟着相应的行块在各处理器之间传送。处理器之间的数据块交换存在如下规律：0号处理器前一个行块(简称前数据块，后一个行块简称后数据块)始终保持不变，将后数据块发送给1号处理器，作为1号处理器的前数据块。同时

15、接收1号处理器发送的后数据块作为自己的后数据块。p-1号处理器首先将后数据块发送给编号为p-2的处理器，作为p-2号处理器的后数据块。然后将前数据块移至后数据块的位置上，最后接收p-2号处理器发送的前数据块作为自己的前数据块。编号为my_rank的其余处理器将前数据块发送给编号为my_rank+1的处理器，作为my_rank+1号处理器的前数据块。将后数据块发送给编号为my_rank-1的处理器，作为my_rank-1号处理器的后数据块。为了避免在通信过程中发生死锁，奇数号处理器和偶数号处理器的收发顺序被错开。使偶数号处理器先发送后接收；而奇数号处理器先将数据保存在缓冲区中，然后接收偶数号处理

16、器发送的数据，最后再将缓冲区中的数据发送给偶数号处理器。数据块传送的具体过程描述如下：算法21.4雅可比法求对称矩阵特征值的并行过程中处理器之间的数据块交换算法输入：矩阵A的行数据块和向量b的数据块分布于各个处理器中输出：在处理器阵列中传送后的矩阵A的行数据块和向量b的数据块Begin对所有处理器my_rank(my_rank=0,，p-1)同时执行如下的算法：/*矩阵A和向量b为要传送的数据块*/(1) if(my-rank=0)then/*0号处理器*/(1.(1) 数据块发送给1号处理器(1.(2) 1号处理器发送来的后数据块作为自己的后数据块endif(2) if(my-rank=p-

17、1)and(my-rankmod2丰0then/*处理器p-1且其为奇数*/(2.1)fori=m/2tom-1do/*将后数据块保存在缓冲区buffer中*/forj=0ton-1dobufferi-m/2,j=ai,jendforendfor(2.(2) fori=m/2tom-1dobufi-m/2=biendfor(2.(3) fori=0tom/2-1do/*将前数据块移至后数据块的位置上*/forj=0ton-1doai+m/2,j=ai,jendforendfor(2.(4) fori=0tom/2-1dobi+m/2=biendfor(2.(5) p-2号处理器发送的数据块作为

18、自己的前数据块(2.(6) uffer中的后数据块发送给编号为p-2的处理器endifif(my-rank=p-1)and(my-rankmod2=0)then/*处理器p-1且其为偶数*/(3.(1) 数据块发送给编号为p-2的处理器(3.(2) fori=0tom/2-1do/*将前数据块移至后数据块的位置上*/forj=0ton-1doai+m/2,j=ai,jendforendfor(3.(3) fori=0tom/2-1dobi+m/2=biendfor(3.(4) p-2号处理器发送的数据块作为自己的前数据块endif(4) if(my-rankwp-1)and(my-rank丰0

19、then/*其它的处理器*/(4.(1) if(my-rankmod2=0)then/*偶数号处理器*/(i)将前数据块发送给编号为my_rank+1的处理器(ii)将后数据块发送给编号为my_rank-1的处理器(ii)接收编号为my_rank-1的处理器发送的数据块作为自己的前数据块(iv)接收编号为my_rank+1的处理器发送的数据块作为自己的后数据块else/*奇数号处理器*/(v)fori=0tom-1do/*将前后数据块保存在缓冲区buffer中*/forj=0ton-1dobufferi,j=ai,jendforendfor(vi)fori=0tom-1dobufi=biend

20、for(vii)接收编号为my_rank-1的处理器发送的数据块作为自己的前数据块Endendif(viii)接收编号为my_rank+1的处理器发送的数据块作为自己的后数据块(ix)将存于buffer中的前数据块发送给编号为my_rank+1的处理器(x)将存于buffer中的后数据块发送给编号为my_rank-1的处理器endif各处理器并行地对其局部存储器中的非主对角元素aj进行消去，首先计算旋转参数并对第i行和第j行两行元素进行旋转行变换。然后通过扩展收集操作将相应的旋转参数及第i列和第j列的列号按处理器编号连接起来并广播给所有处理器。各处理器在收到这些旋转参数和列号之后，按0,1，,

21、p-1的顺序依次读取旋转参数及列号并对其m行中白第i列和第j列元素进行旋转列变换。经过一轮计算的2p-1次数据交换之后，原矩阵A的所有非主对角元素都被消去一次。此时，各处理器求其局部存储器中的非主对角元素的最大元localmax,然后通过归约操作的求最大值运算求得将整个n阶矩阵非主对角元素的最大元max,并广播给所有处理器以决定是否进行下一轮迭代。具体算法框架描述如下：算法21.5雅可比法求对称矩阵特征值的并行算法输入：矩阵Anxn,输出：矩阵A的主对角元素即为A的特征值Begin对所有处理器my_rank(my_rank=0,，p-1)同时执行如下的算法：(a)fori=0tom-1dobi

22、=myrank*m+i/*b记录处理器中的行块数据在原矩阵A中的实际行号*/endfor(b)while(|max|)edo/*max为A中所有非对角元最大的绝对值*/(1)fori=my_rank*mto(my_rank+1)*m-2do/*对本处理器内部所有行两两配对进行旋转变换*/forj=i+1to(my_rank+1)*m-1do(1.(1) imodm,t=jmodm/*i,j为进行旋转变换行(称为主行)的实际行号,r,t为它们在块内的相对行号*/(1.(2) (ar,jw0)then/*对四个主元素的旋转变换*/(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/*对两个主行其余元素的旋转变换*/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/*对两个主列在本处理器内的其余元素的旋转变换*/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保存本处理器主行的行号和旋

25、转参数*/temp10=sin1,temp11=cos1,temp12=(float)i,temp13=(float)jelse(vii)Compute:temp10=0,temp11=0,temp12=0,temp13=0endif(1.(3) 有处理器temp1中的旋转参数及主行的行号按处理器编号连接起来并广播给所有处理器，存于temp2中(1.(4) rrent=0(1.(5) forv=1topdo/*根据temp2中的其它处理器的旋转参数及主行的行号对相关的列在本处理器的部分进行旋转变换*/(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/*进彳T2p-2次处理器间的数据交换，并对交换后处理器中所有行两两配对进

27、行旋车t变换*/(2.(1) ta_exchange()/*处理器间的数据交换*/(2.(2) fori=0tom/2-1doforj=m/2tom-1do(1) if(ai,bjw0)en/*对四个主元素的旋转变换*/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/*对两个主行其余元素的旋转变换*/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/*对本处理器内两个主列的其余元素旋转变换*/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保存本处理器主行的行号和旋转参数*/temp10=sini,temp11=cosi,temp12=(float)bi,temp13=(float)bjendforelse Compute:temp10=0,temp11=0,temp12=0,temp13=0endif(ii)将所有处理器temp1中的旋转参数及主行的行号按处理器编号连接起来并广

30、播给所有处理器，存于temp2中(iii) current=0(iv) forv=1topdo/*根据temp2中的其它处理器的旋转参数及主行的行号对相关的列在本处理器的部分进行旋转变换*/ 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()/*进行一轮中的最后一次处理器间的数据交换，使数据回到原来的位置*/(4)localmax=max/*localmax为本处理器中非对角元最大的绝对值*/(5)fori=0tom-1doforj=0ton-1doif(m*my-rank+i)矜thenEndif(ai,jlocalmax)thenlocalmax=ai,jendifendifendforendfor(6)通过Allreduce操作求出所有处理器中loca

32、lmax的最大值为新的max值endwhile在上述算法中，每个处理器在一轮迭代中要处理2p个行块对，由于每一个行块对含有m/2行，因而对每一个行块对的处理要有(m/2)2=m2/4个行的配对，即消去m2/4个非主对角元素.每消去一个非主对角元素要对同行n个元素和同列n个元素进行旋转变换.由于按行划分数据块，对同行n个元素进行旋转变换的过程是在本处理器中进行的.而对同列n个元素进行旋转变换的过程则分布在所有处理器中进行.但由于所有处理器同时在为其它处理器的元素消去过程进行列的旋转变换，对每个处理器而言，对列元素进行旋转变换的处理总量仍然是n个元素。因此，一个处理器每消去一个非主对角元素共要对2

33、n个元素进行旋转变换。而对一个元素进行旋转变换需要2次乘法和1次加法，若取一次乘法运算时间或一次加法运算时间为一个单位时间，则其需要3个单位运算时间，即一轮迭代的计算时间为Ti=3X2pX2nxm2/4=3n3/p；在每轮迭代中，各个处理器之间以点对点的通信方式相互错开交换数据2p-1次，通信量为mn+m,扩展收集操作n(n-1)/(2p)次，通信量为4,另外有归约操作1次，通信量为1，从而不难得出上述算法求解过程中的总通信时间为：T2=4ts+m(n+1)tw(4p-2)+n(n-1)/p+2t3p-1)+2n(n-1)/p+1tw(p-1)(因此雅可比算法求对矩阵特征值的一轮并行计算时间为

34、Tp=T1+T2。我们的大量实验结果说明，对于n阶方阵，用上述算法进行并行计算，一般需要不超过O(logn)轮就可以收敛。MPI源程序请参见章末附录。1.3求对称矩阵特征值的单侧旋转法1.3.1 单侧旋转法的算法描述求解对称方阵特征值及特征向量的雅可比算法在每次消去计算前都要对非主对角元素选择最大元，导致非实际计算开销很大。在消去计算时，必须对方阵同时进行行、列旋转变换，这称之为双侧旋转(Two-sideRotation)变换。在双侧旋转变换中，方阵的行、列方向都有数据相关关系，使得整个计算中的相关关系复杂，特别是在对高阶矩阵进行特征值求解时，增加了处理器间的通信开销。针对传统雅可比算法的缺点

35、，可以将双侧旋转改为单侧旋转，得出一种求对称矩阵特征值的快速算法。这一算法对矩阵仅实施列变换，使得数据相关关系仅在同列之间，因此方便数据划分，适合并行计算，称为单侧旋转法(One-sideRotation)若A为一对称方阵，入是A的特征值，x是A的特征向量，则有：Ax=及。又A=AT,所以ATAx=Ax=Xx,这说明了是ATA的特征值，x是ATA的特征向量，即ATA的特征值是A的特征值的平方，并且它们的特征向量相同。我们使用18.7.1节中所介绍的Givens旋转变换对A进行一系列的列变换，得到方阵Q使其各列两两正交，即AV=Q这里V为表示正交化过程的变换方阵。因QtQ=(AV)TAV=VTA

36、TAV=diag(B1,B2,，Sn)为n阶对角方阵，可见这里B1,62,，Bn是矩阵ATA的特征值，VI各列是Ata的特征向量。由此可得：6i=&(i=1,2,n)。设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号处理器从其余各处理器中接收子矩阵a和e,得到经过变换的矩阵A和I:/*0号处理器负责计算特征值*/(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分解求矩阵特征值算法的一轮时间复杂度为QR分解、矩阵转置和矩阵相乘算法的时间复杂度之和。1.4.2 QR方法求一般矩阵全部特征值的并行算法并行QR分解求矩阵特征值的思想就是反复运用并行QR分解和并行矩阵相乘算法进行迭代，直到矩阵序列a收敛于一个上三角矩阵或块上三角矩阵为止。具体的并行算法描述如下：算法21.9QR方法求一般矩阵全部特征值的并行算法输入：矩阵An。,单位矩阵Q,输出：矩阵特征值日genvalueBegin对所有处理器my_rank(my_rank=0,，p-1)同时执行如下的算法：(a)while(p)and(c

41、ount0)and(my_rankp)thenp=ai,jendifendforendforendwhile(b)fori=1tondoEigenvalue。=ai,iendfor在上述算法的一轮迭代中，实际上是使用算法18.12、18.1、18.6并行地进行矩阵的QR分解、矩阵的转置、矩阵的相乘各一次，因此，一轮迭代的计算时间为上述算法的时间复杂度之和。MPI源程序请参见所附光盘。1.5小结矩阵的特征值与特征向量在工程技术上应用十分广泛，在诸如系统稳定性问题和数字信号处理的K-L变换中，它都是最基本、常用的算法之一。本章主要讨论了矩阵特征值的几种基本算法，关于该方面的其它讲解与分析读者可参考文献1、2。参考文献1 .陈国良编著.并行算法的设计与分析（修订版）.高等教育出版社，2002.112 .陈国良，陈峻编著.VLSI计算理论与并行算法.中国科学技术大学出版社，1991附录求对称矩阵特征值的雅可比并行算法MPI源程序源程序cjacobi.c#includestdio.h#includestdlib.h#includempi.h#includemath.h#includestring.h#defineE0.00