并行计算矩阵特征值计算Word文档下载推荐.docx

并行计算矩阵特征值计算Word文档下载推荐.docx

《并行计算矩阵特征值计算Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《并行计算矩阵特征值计算Word文档下载推荐.docx（47页珍藏版）》请在冰豆网上搜索。

特征值λi对应的特征向量为xi。

乘幂法的做法是：

①取n维非零向量v0作为初始向量；

②对于k=1,2,„,做如下迭代：

直至uuk=Avk-1vk=uk/║uk║∞

k+1∞

uk∞<

ε为止，这时vk+1就是A的绝对值最大的特征值λ1所对应的特征向

量x1。

若vk-1与vk的各个分量同号且成比例，则λ1=║uk║∞；

若vk-1与vk的各个分量异号且成比例，则λ1=-║uk║∞。

若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时间，则因为一轮计算要做一次矩阵向量相乘、一次求最大元操作和一次规格化操作，所以下述乘幂法串行算法21.1的一轮计算时间为n2+2n=O（n2）。

算法21.1单处理器上乘幂法求解矩阵最大特征值的算法

输入：

系数矩阵An×

n，初始向量vn×

1，ε

输出：

最大的特征值max

Begin

while（│diff│>

ε）do

（1）fori=1tondo

（1.1）sum=0

（1.2）forj=1tondo

sum=sum+a[i,j]*x[j]

endfor

（1.3）b[i]=sum

（2）max=│b[1]│

（3）fori=2tondo

if（│b[i]│>

max）thenmax=│b[i]│endif

（4）fori=1tondo

x[i]=b[i]/max

（5）diff=max-oldmax,oldmax=max

endwhile

End

乘幂法求矩阵特征值由反复进行矩阵向量相乘来实现，因而可以采用矩阵向量相乘的数据划分方法。

设处理器个数为p，对矩阵A按行划分为p块，每块含有连续的m行向量，这里m=⎡n/p⎤，编号为i的处理器含有A的第im至第（i+1）m-1行数据，（i=0,1,„,p-1），初始向量v被广播给所有处理器。

各处理器并行地对存于局部存储器中a的行块和向量v做乘积操作，并求其m维结果向量中的最大值localmax，然后通过归约操作的求最大值运算得到整个n维结果向量中的最大值max并广播给所有处理器，各处理器利用max进行规格化操作。

最后通过扩展收集操作将各处理器中的m维结果向量按处理器编号连接起来并广播给所有处理器，以进行下一次迭代计算，直至收敛。

具体算法框架描述如下：

算法21.2乘幂法求解矩阵最大特征值的并行算法

对所有处理器my_rank（my_rank=0,„,p-1）同时执行如下的算法:

ε）do/*diff为特征向量的各个分量新旧值之差的最大值

*/

（1）fori=0tom-1do/*对所存的行计算特征向量的相应分量*/

（1.2）forj=0ton-1do

（2）localmax=│b[0]│/*对所计算的特征向量的相应分量

求新旧值之差的最大值localmax*/

（3）fori=1tom-1do

localmax）thenlocalmax=│b[i]│endif

（4）用Allreduce操作求出所有处理器中localmax值的最大值max

并广播到所有处理器中

（5）fori=0tom-1do/*对所计算的特征向量归一化*/

b[i]=b[i]/max

（6）用Allgather操作将各个处理器中计算出的特征向量的分量的新值组合并广播到所有处理器中

（7）diff=max-oldmax,oldmax=max

endwhileEnd

若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时

间，一轮迭代的计算时间为mn+2m；

一轮迭代中，各处理器做一次归约操作，通信量为1，一次扩展收集操作，通信量为m，则通信时间为4ts（

（

p-1）+（m+1）tw（p-1）。

因此乘幂法的一轮并行计算时间为Tp=mn+2m+4t-1）+（m+1）tw（p-1）。

MPI源程序请参见所附光盘。

1.2求对称矩阵特征值的雅可比法

若矩阵A=[aij]是n阶实对称矩阵，则存在一个正交矩阵U，使得

UTAU=D

其中D是一个对角矩阵，即

⎡λ1⎢0D=⎢⎢⎢⎢⎣000⎤0⎥⎥⎥⎥λn⎥⎦λ2

这里λi（i=1,2,„,n）是A的特征值，U的第i列是与λi对应的特征向量。

雅可比算法主要是通过正交相似变换将一个实对称矩阵对角化，从而求出该矩阵的全部特征值和对应的特征向量。

因此可以用一系列的初等正交变换逐步消去A的非对角线元素，从而使矩阵A对角化。

设初等正交矩阵为R（p,q,θ），其（p,p）（q,q）位置的数据是cosθ，（p,q）（q,p）位置的数据分别是-sinθ和sinθ（p≠q），其它位置的数据和一个同阶数的单位矩阵相同。

显然可以得到：

R（p,q,θ）TR（p,q,θ）=In

不妨记B=R（p,q,θ）TAR（p,q,θ），如果将右端展开，则可知矩阵B的元素与矩阵A的元素之间有如下关系：

bpp=appcos2θ+aqqsin2θ+apqsin2θbqq=appsin2θ+aqqcos2θ-apqsin2θ

bpq=（（aqq-app）sin2θ）/2+apqcos2θbqp=bpq

bpj=apjcosθ+aqjsinθbqj=-apjsinθ+aqjcosθ

bip=aipcosθ+aiqsinθbiq=-aipsinθ+aiqcosθ

bij=aij

其中1≤i,j≤n且i,j≠p,q。

因为A为对称矩阵，R为正交矩阵，所以B亦为对称矩阵。

若要求矩阵B的元素bpq=0，则只需令（（aqq-app）sin2θ）/2+apqcos2θ=0，即：

tg2θ=

-apq（aqq-app2

在实际应用时，考虑到并不需要解出θ，而只需要求出sin2θ，sinθ和cosθ就可以了。

如果限制θ值在-π/2<

2θ≤π/2，则可令

m=-an=（aqq-apq,pp）,2容易推出：

wsin2θ=w,sinθ=,2（1+-w2）1w=sgn（n）mm2+n2cosθ=-sin2θ利用sin2θ，sinθ和cosθ的值，即得矩阵B的各元素。

矩阵A经过旋转变换，选定的非主

对角元素apq及aqp（一般是绝对值最大的）就被消去，且其主对角元素的平方和增加了2a2pq，而非主对角元素的平方和减少了2a2pq，矩阵元素总的平方和不变。

通过反复选取主元素apq，并作旋转变换，就逐步将矩阵A变为对角矩阵。

在对称矩阵中共有（n2-n）/2个非主对角元素要被消去,而每消去一个非主对角元素需要对2n个元素进行旋转变换,对一个元素进行旋转变换需要2次乘法和1次加法，若各取一次乘法运算时间或一次加法运算时间为一个单位时间，则消去一个非主对角元素需要3个单位运算时间，所以下述算法21.3的一轮计算时间复杂度为（n2-n）/2*2n*3=3n2（n-1）=O（n3）。

算法21.3单处理器上雅可比法求对称矩阵特征值的算法

矩阵An×

n，ε

矩阵特征值Eigenvalue

（1）while（max>

（1.1）max=a[1,2]

（1.2）fori=1tondo

forj=i+1tondo

if（│a[i,j]）│>

max）thenmax=│a[i,j]│，p=i，q=jendif

（1.3）Compute:

m=-a[p,q]，n=（a[q,q]-a[p,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），

b[p,p]=a[p,p]*co1*co1+a[q,q]*si1*si1+a[p,q]*si2，

b[q,q]=a[p,p]*si1*si1+a[q,q]*co1*co1-a[p,q]*si2，

b[q,p]=0，b[p,q]=0

（1.4）forj=1tondo

if（（j≠p）and（j≠q））then（i）b[p,j]=

a[p,j]*co1+a[q,j]*si1（ii）b[q,j]=

-a[p,j]*si1+a[q,j]*co1

endif

（1.5）fori=1tondo

if（（i≠p）and（i≠q））then（i）b[i,p]=a[i,p]*co1+a[i,q]*si1（ii）b[i,q]=-a[i,p]*si1+a[i,q]*co1endifendfor（1.6）fori=1tondoforj=1tondoa[i,j]=b[i,j]endforendforendwhile

（2）fori=1tondoEigenvalue[i]=a[i,i]endforEnd

串行雅可比算法逐次寻找非主对角元绝对值最大的元素的方法并不适合于并行计算。

因此，在并行处理中，我们每次并不寻找绝对值最大的非主对角元消去，而是按一定的顺序将A中的所有上三角元素全部消去一遍，这样的过程称为一轮。

由于对称性，在一轮中，A的下三角元素也同时被消去一遍。

经过若干轮，可使A的非主对角元的绝对值减少，收敛于一个对角方阵。

具体算法如下：

设处理器个数为p，对矩阵A按行划分为