《矩阵特征值计算.pdf》由会员分享，可在线阅读，更多相关《矩阵特征值计算.pdf（32页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、9矩阵特征值计算在实际的工程计算中，经常会遇到求n 阶 方 阵 A 的特征值(Ei g e n v a l u e)与特征向量(Ei g e n v e c t o r)的问题。对于-,个方阵/,如果数值A,使方程组Ax=Xx即(4力”)尸0有非零解向量(S o l u t i o n V e c t o r ,则称力为方阵/的特征值，而非零向量x为特征值义所对应的特征向量，其 中I”为阶单位矩阵。由于根据定义直接求矩阵特征值的过程比较复杂，因此在实际计算中，往往采取一些数值方法。本章主要介绍求一般方阵绝对值最大的特征值的乘幕(P o w e r)法、求对称方阵特征值的雅可比法和单侧旋转(On
2、 e-si d e R o t a t i o n)法以及求一般矩阵全部特征值的Q R方法及一些相关的并行算法。1.1 求解矩阵最大特征值的乘幕法1.1.1 乘塞法及其串行算法在许多实际问题中，只需要计算绝对值最大的特征值,而并不需要求矩阵的全部特征值。乘幕法是一种求矩阵绝对值最大的特征值的方法。记 实 方 阵 力 的 个 特 征 值 为 4z=(l,2,且满足：|2,|22|23|-|2|特征值不对应的特征向量为X*乘基法的做法是：取维非零向量V。作为初始向量；对于Q1,2,，做如下迭代：uk=Avk.%=殊/|”48直至卜|以+1 1 M L s)do(l)for z=l to n do(
3、l.l)5W /W=0(1.2)fo r j=1 to dosum=sum+aij*x/end for(1.3)/z=sumend for(2)max=h|(3)for i=2 to n doif(|bi|max)then max=|bi|end ifend for(4)for i=1 to n doxz=bihnaxend for(5)di=tnax-oldmax,oldmax=niaxend whileEnd1.1.2乘鬲法并行算法乘募法求矩阵特征值由反复进行矩阵向量相乘来实现，因而可以采用矩阵向量相乘的数据划分方法。设处理器个数为p,对矩阵A按行划分为p块，每块含有连续的m行向量，这里,
4、编号为i 的处理器含有X 的第加至第行数据，(/=0,1,p-1),初始向量v 被广播给所有处理器。各处理器并行地对存于局部存储器中0 的行块和向量v 做乘积操作，并求其m维结果向量中的最大值/。想加办，然后通过归约操作的求最大值运算得到整个维结果向量中的最大值加水并广播给所有处理器，各处理器利用，以进行规格化操作。最后通过扩展收集操作将各处理器中的加维结果向量按处理器编号连接起来并广播给所有处理器，以进行下一次迭代计算，直至收敛。具体算法框架描述如下：算法2 1.2 乘幕法求解矩阵最大特征值的并行算法输入：系数矩阵4 x“，初始向量V”xi，E输出：最大的特征值7以Begin对所有处理器my
5、_rank(my_rank=O,pl)同时执行如下的算法：while(|diff|d。产d iff为特征向量的各个分量新旧值之差的最大值*/(l)for Z=0 to m-do 对所存的行计算特征向量的相应分量*/(1.1)sum=0(1.2)fo r/=0 to n-dosum=suni+ai2/end for(1.3)Z z=sumend for(2localmax=ft0|/*对所计算的特征向量的相应分量求新旧值之差的最大值localmax*/(3)for z=l to m-doif(|bi|localmax)then locahnax=bi end ifend for(4)用 Allr
6、educe操作求出所有处理器中locabnax值的最大值tnax并广播到所有处理器中(5)for z=Oto m-do/*对所计算的特征向量归一化*/hi=bi/maxend for(6)用 Allgather操作将各个处理器中计算出的特征向量的分量的新值组合并广播到所有处理器中(7)diff=max-oldmax,oldmax=maxend whileEnd若各取一次乘法和加法运算时间、次除法运算时间、一次比较运算时间为一个单位时间，一轮迭代的计算时间为m+2如一轮迭代中，各处理器做一次归约操作，通信量为1,一次扩展收集操作，通信量为相，则通信时间为4G(V 7-D +(机+1)9 5-1)
7、。因此乘塞法的一轮并行计算时间为 7P=mn+2m+4ts(yp-1)+(w+l)/,v(p-1)。MPI源程序请参见所附光盘。1.2求对称矩阵特征值的雅可比法1.2.1 雅可比法求对称矩阵特征值的串行算法若矩阵/=%是阶实对称矩阵，则存在一个正交矩阵u,使得UTAU=D其中。是一个对角矩阵，即40A0D=0丸 2A0MM 0M00A几 这里小(i=l,2,是4 的特征值，。的第i 列是与不对应的特征向量。雅可比算法主要是通过正交相似变换将一个实对称矩阵对角化，从而求出该矩阵的全部特征值和对应的特征向量。因此可以用 系列的初等正交变换逐步消去/的非对角线元素，从而使矩阵N对角化。设初等正交矩阵
8、为R(p,q,0),其(p,p)(q,q)位置的数据是cos/(p,q)(q,p)位置的数据分别是-sinJ和 sin O S,g),其它位置的数据和一个同阶数的单位矩阵相同。显然可以得到：R(P，q6)TR(P，q,0)=In不妨记B=R(p,q,9)TAR(p,q,6),如果将右端展开，则可知矩阵B 的元素与矩阵A 的元素之间有如下关系：bpp=。川 威 升%/油 2火他gsin29bpq=（3 卯-%）sin2（9）/2+a的 cos26bPj=a0 cos 什&7sin。b cos火arsinebu=aiibqq=appsm23+aqq cos23-aP9sin20bqp=愤bqj=-
9、aPjSinO+a c os8biq=-a3sin。+acos。其 中 i s i j s 且 因 为/为 对 称 矩 阵，R 为正交矩阵，所以8 亦为对称矩阵。若要求矩阵B 的元素bpq=0,则只需令(aw-app)sm26)l2+apil cos20=O,即：(aqq-Opp)2在实际应用时，考虑到并不需要解出仇而只需要求出s i n 2，s i n。和c o s。就可以了。如果限制9值在皿Z 2 V 2 0 WJ T/2,则可令加=-“掰,=皿),.=s g n(n)I 丁 ,27 m +n容易推出：WI；2s i n 2 =w,s i n 0=.,c o s。=M l-s i n 0J
10、 2(l +J l _ w2)利 用s i n 2 as i n。和c o s d的值，即得矩阵8的各元素。矩阵/经过旋转变换，选定的非主对角元素a及。伊(一般是绝对值最大的)就被消去，且其主对角元素的平方和增加了2%,而非主对角元素的平方和减少了 2/%，矩阵元素总的平方和不变。通过反复选取主元 素 并 作 旋 转 变 换，就逐步将矩阵Z变为对角矩阵。在对称矩阵中共有(2-)/2个非主对角元素要被消去，而每消去一个非主对角元素需要对2 个元素进行旋转变换，对一个元素进行旋转变换需要2次乘法和1次加法，若各取一次乘法运算时间或一次加法运算时间为一个单位时间，则消去一个非主对角元素需要3个单位运
11、算时间，所以下述算法2 1.3的轮计算时间复杂度为(2-)/2*2 *3=3 2(-1尸0(3)。算法2 1.3 单处理器上雅可比法求对称矩阵特征值的算法输入：矩阵N”X ，输出：矩阵特征值EigenvalueBegin(l)while(m ax 8)do(1.1)/wt z x=a l,2(1.2)fo r i=to n dofor j=z+1 to n doif(|max)then max=|aij ,p=i,q=j end ifend forend for(1.3)C o i n pu t e:m=-ap.q,n=(aq,q-ap,p/2,片s g n()*优qr t(例*机+*)，s
12、i 2=w,s=w/s qr t(2(l+s qr t(l-w*w),c o l=s qr t(l-5 z l*5 z l),bp,p=/,q=-叩,p*s+叩,q*colend ifend for(1.6)for i=l to n dofor j=l to n doai,j=bi,jend forend forend while(2)for z=l to n doEigenvaluei=ai,zend forEnd1.2.2雅可比法求对称矩阵特征值的并行算法串行雅可比算法逐次寻找非主对角元绝对值最大的元素的方法并不适合于并行计算。因此，在并行处理中，我们每次并不寻找绝对值最大的非主对角元消去
13、，而是按一定的顺序将力中的所有上三角元素全部消去一遍，这样的过程称为一轮。由于对称性，在一轮中，力的下三角元素也同时被消去一 遍。经过若干轮，可 使A的非主对角元的绝对值减少，收敛于一个对角方阵。具体算法如下：设处理器个数为p,对 矩 阵A按行划分为2 p块，每块含有连续的m/2行向量，记m=nlp,这 些 行 块 依 次 记 为 川，并将山,与4存放与标号为，的处理器中。每轮计算开始，各处理器首先对其局部存储器中所存的两个行块的所有行两两配对进行旋转变换，消去相应的非对角线元素。然后按图2 L 1 所示规律将数据块在不同处理器之间传送，以消去其它非主对角元素。图 1.1 p=4时的雅可比算法
14、求对称矩阵特征值的数据交换图这里长方形框中两个方格内的整数被看成是所移动的行块的编号。在要构成新的行块配对时，只要将每个处理器所对应的行块按箭头方向移至相邻的处理器即可，这样的传送可以在行块之间实现完全配对。当编号为i 和/的两个行块被送至同处理器时，令编号为i 的行块中的每行元素和编号为J 的行块中的每行元素配对，以消去相应的非主对角元素，这样在所有的行块都两两配对之后，可以将所有的非主对角元素消去一遍，从而完成一轮计算。由 图 1.1可以看出，在一轮计算中，处理器之间要2 p-l次交换行块。为保证不同行块配对时可以将原矩阵A的非主对角元素消去，引入变量b记录每个处理器中的行块数据在原矩阵”
15、中的实际行号。在行块交换时，变量6 也跟着相应的行块在各处理器之间传送。处理器之间的数据块交换存在如下规律：0 号处理器前一个行块(简称前数据块，后个行块简称后数据块)始终保持不变，将后数据块发送给1 号处理器，作 为 1 号处理器的前数据块。同时接收1 号处理器发送的后数据块作为自己的后数据块。p-号处理器首先将后数据块发送给编号为p-2的处理器，作为p-2号处理器的后数据块。然后将前数据块移至后数据块的位置上，最后接收p-2号处理器发送的前数据块作为自己的前数据块。编 号 为 my_rank的其余处理器将前数据块发送给编号为my_rank+l的处理器，作为my_rank+l号处理器的前数据
16、块。将后数据块发送给编号为my_rank-l的处理器，作为my_rank-l号处理器的后数据块.为了避免在通信过程中发生死锁，奇数号处理器和偶数号处理器的收发顺序被错开。使偶数号处理器先发送后接收；而奇数号处理器先将数据保存在缓冲区中，然后接收偶数号处理器发送的数据，最后再将缓冲区中的数据发送给偶数号处理器。数据块传送的具体过程描述如下：算法2 1.4 雅可比法求对称矩阵特征值的并行过程中处理器之间的数据块交换算法输入：矩阵Z 的行数据块和向量6 的数据块分布于各个处理器中输出：在处理器阵列中传送后的矩阵N的行数据块和向量b的数据块Begin对所有处理器my_rank(my_rank=O,p-
17、1)同时执行如下的算法：/*矩阵A和向量b为要传送的数据块*/(1)if(my-rank=O)then/*0 号处理器*/(1.1)将后 数据块发送给1号处理器(L2)接 收 1 号处理器发送来的后数据块作为自己的后数据块end if(2)if(my-rank=p-1)and(my-rank mod 2 声 0)then/*处理器p-且其为奇数*/(2.1)for z=w/2 to w-1 do/*将后数据块保存在缓冲区buffer中*/for j=0 to n-1 dobufferi-m/2j=aijend forend for(2.2)fo r i=mll to m-dohufi-ni/2
18、=biend for(2.3)fo r i=0 to m/2-1 do/*将前数据块移至后数据块的位置上*/for j=0 to n-1 doai+m/2j=aijend forend for(2.4)fo r z=0 to w/2-1 dobi+m/2=hiend for(2.5)接收 p-2号处理器发送的数据块作为自己的前数据块(2.6)将b uffer中的后数据块发送给编号为p-2的处理器end if(3)if(my-rank=/?-l)and(my-rank mod 2=0)then/*处理器 pl 且其为偶数*/(3.1)将后 数据块发送给编号为p-2 的处理器(3.2)fo r z
19、=0 to zw/2-1 do/*将前数据块移至后数据块的位置上*/fory=0 to n-doai+m/2j=aijend forend for(3.3)fo r z=0 to m/2-1 dobi+m/2=biend for(3.4)接收 p-2号处理器发送的数据块作为自己的前数据块end if(4)if(my-rank，and(my-rank*0)then/*其它的处理器*/(4.1)if(my-rank mod 2=0)then/*偶数号处理器*/将前数据块发送给编号为myank+l的处理器(ii)将后数据块发送给编号为my_rank-l的处理器(ii)接收编号为my_rank-l的处
20、理器发送的数据块作为自己的前数据块(iv)接收编号为my_rank+l的处理器发送的数据块作为自J 的后数据块else/*奇数号处理器*/(v)for i=0 to w-1 do/*将前后数据块保存在缓冲区buffer 11*/for7=0 to n-dobufferij=aijend forend for(vi)for z=0 to m A dobufi=biend for(vii)接收编号为my_rank-l的处理器发送的数据块作为自己的前数据块(viii)接收编号为my_rank+l的处理器发送的数据块作为自己的后数据块(ix)将存于儿娘尸中的前数据块发送给编号为my_rank+l的处理
21、器(x)将存于buffer中的后数据块发送给编号为my_rank-l的处理器end ifend ifEnd各处理器并行地对其局部存储器中的非主对角元素的 进行消去，首先计算旋转参数并对第i 行和第j 行两行元素进行旋转行变换。然后通过扩展收集操作将相应的旋转参数及第/列和第j 列的列号按处理器编号连接起来并广播给所有处理器。各处理器在收到这些旋转参数和列号之后，按 0,1,夕-1 的顺序依次读取旋转参数及列号并对其加行中的第，列和第/列元素进行旋转列变换。经过一轮计算的2/1 次数据交换之后，原矩阵力的所有非主对角元素都被消去一次。此时，各处理器求其局部存储器中的非主对角元素的最大元/。馆加然
22、后通过归约操作的求最大值运算求得将整个阶矩阵非主对角元素的最大元机衣，并广播给所有处理器以决定是否进行下轮迭代。具体算法框架描述如下：算法2 1.5 雅可比法求对称矩阵特征值的并行算法输入：矩阵输出：矩阵力的主对角元素即为力的特征值Begin对所有处理器my_rank(my_rank=O,p-1)同时执行如卜的算法：(a)for z=0 to/w-1 dobi=myrank*/w+z/*b 记录处理器中的行块数据在原矩阵A 中的实际行号*/end for(b)while(|max e)do/*m a x为 A 中所有非对角元最大的绝对值*/(1)for z=my_rank*/to(my_ran
23、k+l)*m-2 do/*对本处理器内部所有行两两配对进行旋转变换*/forj=/+l to(my_rank+l)*”l do(1.1)7 T mod m,t=j mod m/*i,j为进行旋转变换行(称为主行)的实际行号，厂 为它们在块内的相对行号*/(1.2)if(a 力+0)then/*对四个主元素的旋转变换*/(i)Compute:户-J,g=(。犷卜=sg(g)7/(Z*Ag*g),sin2=h,sin 1 =h!sqrt(2(1 +sqrt(1 -7?*/?),cos 1 =sqrt(1 -sin 1 *sin 1),bpp=ar,i*cos 1 *cos 1 *sm 1 +arj
24、*sin2,bqq=arj*s加cos 1 *cos 1 -arj*sin2,bpq=0,bqp=。(ii)for v=0 to/?-l do/*对两个主行其余元素的旋转变换*/if(v*i)and(u/)thenbry=ar,v*cosl+at,v*sinat9v=-ar,v*sin+a/,v*cosend ifend for(iii)for v=0 to n-doif(v*/)and(u#/)thenar9v=brvend ifend for(iv)for v=0 to m-do/*对两个主列在本处理器内的其余元素的旋转变换*/if(v/r)and(v /)thenbiv=avj*cos+
25、ayj*sinl仇 匕/=匕 4*sin+Q 匕/*coslend ifend for(v)for v=0 to/-ldoif(v*r)and(v/)thenav,i=bivend ifend for(vi)Compute:ar,i=bpp,atj=bqq,arj=bpq,at,i=bqp,/*用 tempi保存本处理器主行的行号和旋转参数*/temp 1 0=s加 1,temp I 1 =cos 1,temp 1 2=(float)z,temp 1 3=(float)/else(vii)Compute:temp 1 0=0,temp 11=0/temp 12=0,temp 13=0end i
26、f(1.3)将所有处理器tem pi中的旋转参数及主行的行号按处理器编号连接起来并广播给所有处理器，存于tem pi中(1 A)current=0(1.5)fo r v=l to p do/*根据拒叫2 中的其它处理器的旋转参数及主行的行号对相关的列在本处理器的部分进行旋转变换*/(i)Compute:s 1 =temp2(y-1 )*4+0,c 1 =temp2(y-1 )*4+1,/I=(int)/e/np2(v-1 )*4+2,j 1 =(int)/ezwp2(v-1 )*4+3(ii)if(sl、c l、il、jl 中有一不为 0)thenif(my-rank/current then
27、for z=0 to m-doziz=az,il*cl+azjl*sla zjl=-az,i*s+azg*cend forfor z=0 to/-I doazj=zizend forend ifend if(in)current=current+1end forend forend for(2)for counter=1 to 2P-2 do/*进 行 2p-2次处理器间的数据交换，并对交换后处理器中所有行两两配对进行旋转变换*/(2.1 )Data_exchange()/*处理器间的数据交换*/(2.2)fo r z=0 to mJ2-dofor j=niH to do(i)if(6/,/+
28、0)then/*对四个主元素的旋转变换*/Compute:户-叩,则，g=(如 即 -叩向/2，h=sgn(g)*f/sqt 仁户g*g)，sin2=h,sin I=A/sq片(2*(1 +sqrt(1 /*),cos 1 =sqrt 1 -sin*sin 1),bpp=ai,bicosl*cosl+af,bf*sinl*sinl+ai,bj*sin2fbqq=ai,bi*s加l*s加l+此 /*cos 1 *cos 1 -ai,bj*sin2,bpq=0,bqp=Gfor v=0 to n-do/*对两个主行其余元素的旋转变换*/if(v*hi)and(u 彳 /)thenbrv=ai,v*
29、cos+aj,v*sinaJ,v=-af,v*sin+叫，0*coslend ifend forfor v=0 to n-doif(v丰 口)and(v/)thenaiyv=brvend ifend forfor v=0 to w-1 do/*对本处理器内两个主列的其余元素旋转变换*/if 声 /)and(v 彳 J)thenbiv=av,hi*cos+ay,bj*sinav,b/=-av,/?/*sinl+av,6/*cosend ifend forfor v=0 to m-doif(v/I)and(v 丰j)then力,bi=bivend ifend forCompute:bi=bpp,a
30、f，b/=bqq,。口，hU=hP(l，bi=bqp/*用 Zempl保存本处理器主行的行号和旋转参数*/temp 1 0=5/H 1 ,temp 1 1 =cos 1,temp 1 2=(float)/?z,templ3=(float)/elseCompute:temp 1 0=0,temp 11=0,tempi 2=09temp3=0end if(ii)将所有处理器tem pi中的旋转参数及主行的行号按处理器编号连接起来并广播给所有处理器，存于tem pi中(iii)current=0(iv)for v=l top do/*根据加w 2 中的其它处理器的旋转参数及主行的行号对相关的列在本处
31、理器的部分进行旋转变换*/Compute:s 1 =templ(v-1)*4+0,c 1 =temp2(v-1)*4+1,/I=(int)/ewp2(v-1 )*4+2,jl=(int)/ewp2(v-1 )*4+3 if(si、c l、中有一不为 0)thenif(my-rank*current)thenfor z=0 to m-doziz=az9i*cl+azj*vla zji=-az,ils+azj*clend forfor z=0 to m-doaz,i=zizend forend ifend ifcurrent=current+1end forend forend forend fo
32、r(3)Data-exchange()/*进行一轮中的最后一次处理器间的数据交换，使数据回到原来的位置*/()locahnax=inax localmax为本处理器中非对角元最大的绝对值*/(5)for z=0 to 777-1 dofory=0 to/?-l doif(w*my-rank+f)#j)thenif(|aij|localmax)then localmax=end ifend ifend forend for(6)通过Allreduce操作求出所有处理器中localmax的最大值为新的max值end whileEnd在上述算法中，每个处理器在-轮迭代中要处理2 2 个行块对，由 于
33、 每 个行块对含有m12行，因而对每一个行块对的处理要有(m/2)2=/4个行的配对，即消去苏/4 个非主对角元素.每消去一个非主对角元素要对同行个元素和同列n个元素进行旋转变换.由于按行划分数据块，对 同 行n个元素进行旋转变换的过程是在本处理器中进行的.而对同列个元素进行旋转变换的过程则分布在所有处理器中进行.但由于所有处理器同时在为其它处理器的元素消去过程进行列的旋转变换，对每个处理器而言，对列元素进行旋转变换的处理总量仍然是个元素。因此，一 个处理器每消去一个非主对角元素共要对2个元素进行旋转变换。而对一个元素进行旋转变换需要2 次乘法和1 次加法，若取一次乘法运算时间或一次加法运算时
34、间为一个单位时间，则其需要3 个单位运算时间，即一轮迭代的计算时间为方=3X 2p X2n 乂舟4=3曲p；在每轮迭代中，各个处理器之间以点对点的通信方式相互错开交换数据2p-l次，通信量为机+如 扩展收集操作(-l)/(2p)次，通信量为4,另外有 归 约 操 作 1 次，通 信 量 为 1,从而不难得出上述算法求解过程中的总通信时间为：T2=%+m(n +T)tw(4p-2)+(-1)/p +2 k (后 一 1)+2 n(-l)/p +ltw(p 一 D因此雅可比算法求对矩阵特征值的一轮并行计算时间为7；=q +r2o 我们的大量实验结果说明，对于”阶方阵，用上述算法进行并行计算，一般需
35、要不超过。(/og)轮就可以收敛。MPI源程序请参见章末附录。1.3求对称矩阵特征值的单侧旋转法1.3.1 单侧旋转法的算法描述求解对称方阵特征值及特征向量的雅可比算法在每次消去计算前都要对非主对角元素选择最大元，导致非实际计算开销很大。在消去计算时，必须对方阵同时进行行、列旋转变换，这称之为双侧旋转(Two-side Rotation)变换。在双侧旋转变换中，方阵的行、列方向都有数据相关关系，使得整个计算中的相关关系复杂，特别是在对高阶矩阵进行特征值求解时，增加了处理器间的通信开销。针对传统雅可比算法的缺点，可以将双侧旋转改为单侧旋转，得 出 一 种求对称矩阵特征值的快速算法。这一 算法对矩
36、阵仅实施列变换，使得数据相关关系仅在同列之间，因此方便数据划分，适合并行计算，称为单侧旋转法(One-side Rotation)。若/为一对称方阵，大是/的特征值，x 是/的特征向量，则有：心=&。5LA=AT,所以ArA x=A U=U x,这说 明 了 好 是/，的特征值，x 是 的 特 征 向 量，即的特征值是力的特征值的平方，并且它们的特征向量相同。我们使用18.7.1节中所介绍的Givens旋转变换对A 进行一系列的列变换，得到方阵。使其各列两两正交，即这 里/为表示正交化过程的变换方阵。因。/。产=diog(%,力，,6“)为阶对角方阵，可见这里久,力，力是矩阵工，的特征值，忆的
37、各列是“，的特征向量。由此可得：0=4 2 (i=l,2,)。设。的第，列 为%，则当，|*|2=石K=亚=同。因此在将A 进行列变换得到各列两两正交的方阵。后，其各列的谱范数|明12即为力的特征值的绝对值。记 H的第，列为匕，A V=Q,所以Avi=qi。又因为修为力的特征向量，所以力为=无修，即推得无匕二%。因此无的符号可由向量名及均的对应分量是否同号判别，实际上在具体算法中我们只要判断它们的第一个分量是否同号即可。若相同，则取正值，否则取负值。求对称矩阵特征值的单侧旋转法的串行算法如下：算法2 1.6求对称矩阵特征值的单侧旋转法输入：对称矩阵力，精度e输出：矩阵4 的特征值存于向量8 中
38、Begin(l)while(pe)p=0for i=l to n dofor j=i+to n do(1.1)s m0=0,sum 1 =0,sum2=0(1.2)for k=to n dosum0=sum0+ak,z*ak,jsum=suml+akti*ak,i5WW2=sum2+ak,jY 0 人,力end for(1.3)if(|suni0|)then(i)aa=2*sum0bb=suml-suni2rr=sqrt(aa*aabb*bb)(ii)if(hb0)thenc=sqrt(防+AT)/(2*/T)s=aa/(2*c)else5=sqrt(rr-/6)/(2*rr)c=aa/(2*
39、rr*s)end if(iii)for k=to n dotefnpk=c*ak,i+s*ak,jakj=-s*ak,i+c*akjend for(iv)for k=to n dotempk=c*ek,i+s*ekgekj=-s*ek,i+c*ekjend for(v)for k=to n doak,i=tempkend for(vi)for k=to n doek,i=tempkend for(vii)if(|sum0|p)then p=sum0|end ifend ifend forend forend while(2)for z=l to n dosu=0fory=l to n dosu
40、=sua/,i*a/Jend for8/=sqrts if(e0j*a0 j e)then/*各个处理器并行进行对i和/列正交化*/aa=2*s 皿 0励1 -sum2/7=sqrt(m*QQ+bb*bb)if(bb=0)thenc=sqrt(ft/)+rr)/(2*rr)s=aa/(2*rr*c)else5=sqrt(rr-66)/(2*rr)c=aa/(2*rr*s)end iffor v=0 to q-dotempy=cay,i+s*av,jH%/=(-s)*a 匕 i+c*avjend forfor v=0 to z-1 dotemplv=c*ev,is*ev,jend forfor
41、v=0 to q-doav4=tempvend forfor v=0 to z-1 doev,i=tempvend for 仁k+1end ifend forend forend while0 号处理器从其余各处理器中接收子矩阵a 和 e,得到经过变换的矩阵A 和 I：/*0号处理器负责计算特征值*/(b)for i=l to n do(l)s=O(2)for7=1 to n dosz尸sz/+4 /；i *J /,/end for(3)S/=sqrt(5 M)(4)if(IO,j*aO,j 8)a n d (c o，v l 0 0 0)do(1.1 )p=0,count=count+1(1.
42、2)Q R_ D E C O M P O SI T I O N(Z,0,R)/*算法 1 8.1 1 对矩阵 4 进行 Q R 分解*/(1.3)MATRIX_TRANSPOSITION(0/*由于算法 18.11 对 A 进行 QR 分解只 得 到 因 而 要 使 用 算 法 1 8.1 对矩阵。进行转置，得到(2)T=(。午=Q*/(1.4M=MATRIX_PRODUCT(7?,0/*算法 18.5 将矩阵 AQ 相乘，得到新的A矩 阵*/(1.5)for;=1 t o n dofor j=l t o i-1 doif(|aij|p)then p=aij end ifend forend
43、forend while(2)for z=l t o n doEigenvaluei=ai,iend forEnd显然，QR分解求矩阵特征值算法的一轮时间复杂度为QR分解、矩阵转置和矩阵相乘算法的时间复杂度之和。1.4.2 QR方法求一般矩阵全部特征值的并行算法并行QR分解求矩阵特征值的思想就是反复运用并行QR分解和并行矩阵相乘算法进行迭代，直到矩阵序列 4 收敛于一个上三角矩阵或块上三角矩阵为止。具体的并行算法描述如下：算法21.9 QR方法求一般矩阵全部特征值的并行算法输入：矩阵4x“,单位矩阵0，e输出：矩阵特征值EigenvalueBegin对所有处理器m y _ r an k(m y
44、 _ r an k=O,p-1)同时执行如下的算法：(a)while(p e)an d(c o w w/0)an d(m y _ r an k P)then p=aij end ifend forend forend while(b)for z=l to n doEigenvaluei=t/z,zend forEnd在上述算法的一轮迭代中，实际上是使用算法18.12、18.1、18.6并行地进行矩阵的QR分解、矩阵的转置、矩阵的相乘各一次，因此，一轮迭代的计算时间为上述算法的时间复杂度之和。MPI源程序请参见所附光盘。1.5小结矩阵的特征值与特征向量在工程技术上应用十分广泛，在诸如系统稳定性问
45、题和数字信号处理的K-L变换中，它都是最基本、常用的算法之一。本章主要讨论了矩阵特征值的几种基本算法，关于该方面的其它讲解与分析读者可参考文献I、2参考文献 1.陈国良编著.并行算法的设计与分析（修订版）.高等教育出版社，2002.11 2,陈国良，陈 峻 编 著.VLSI计算理论与并行算法.中国科学技术大学出版社，1991附录求对称矩阵特征值的雅可比并行算法MPI源程序源程序cjacobi.c#include*stdio.h#includc stdlib.h#include mpi.h*#includc math.h#include string.hM#dcfincE 0.00001#def
46、ine min-1#dcfinc intsizc sizcoffint)#define charsize sizeof(char)#dcfinc floatsizc sizcofffloat)/*A 是 个 N*N的对称矩阵*/#dcfinc A(x,y)Ax*N+y/*1是一个N*N的单位矩阵*/#definc I(x,y)Ix*N+ydefine a(x,y)ax*N+y/define e(x,y)ex*N+ydefine b(x)bx#dcfinc buflfcr(x,y)buffcrx*N+ydefine buffee(x,y)buffeex*N+yint M,N;int*b;int
47、m,p;int myid,group_size;float*A,*I;float*a,*e;float max;float sum;MPI Status status;float sgn(float v)(float f;if(v=0)f=l;else f=-l;return f;)void rcad_filcA()(int ij;FILE*R1A;3A=fbpen(dataln.lxtTT)；fscanf(fdA,M%d%d”，&M,&N);iRM!=N)(puts(nThe input is error!);cxit(0);)m=N/p;if(N%p!=0)m+;A=(float*)mal
48、loc(floatsize*N*m*p);I=(float*)malloc(floatsizc*N*m*p);fbr(i=0;i M;i-H-)fdr(j=O;jM;j-H-)抬 canf(fdA,”F,A+i*M+j);fclose(fdA);printflInput of file udataIn.txtHnM);printf(M%dt%dnn,M,N);for(i=0;iM;i+)fbr(j=O;jN;j+)printf(%ft,A(i,j);printf(Mnn);for(i=0;iN;i+)for(j=O;jN;j+)if(i=j)Ki,j)=l；else I(i,j)=0;void
49、 send_a()(inti;for(i=l;ip;i+)MPI_Send(&(A(m*i,O),m*N,MPI_FLOAT,i,i,MPI_COMM_WORLD);MPI_Send(&(l(m*i,O),m*N,MPI_FLOAT,i,i,MPI_COMM_WORLD);frcc(A);free(I);void gct_a()(int ij;if(myid=0)(fbr(i=O;im;i-H-)fbr(j=O;jN;j-H-)(a(i,j)=A(ij);e(i,j)=I(i,j);)elseMPI_Recv(a,m*N,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&sta
50、tus);M PI_Recv(e,m*N,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status);int main(int argc,char*argv)float*c;int k;int loop;int ij,v,z,r,t,y；int il,jl;float f,g,h;float sinl,sin2,cosl;float si,cl;float*br,*bt,*bi,*bj,*zi,*zj;float*temp 1 ,*temp2,*buffer,*buflee;int countercurrent;int*buf;intfloat bpp,bqq,bpq,b