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

9 矩阵特征值计算

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

1.1 求解矩阵最大特征值的乘幂法

1.1.1 乘幂法及其串行算法

在许多实际问题中，只需要计算绝对值最大的特征值，而并不需要求矩阵的全部特征值。 乘幂法是一种求矩阵绝对值最大的特征值的方法。记实方阵A的n个特征值为λi i=(1,2, …,n)， 且满足：
│λ1 │≥│λ2 │≥│λ3 │≥…≥│λn │
特征值λi 对应的特征向量为xi 。乘幂法的做法是：①取n维非零向量v0 作为初始向量；②对于k=1,2, …,做如下迭代：

直至 u


k 1 ¥ - uk
uk =Avk-1 vk = uk /║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 ，初始向量v n×1 ，ε
输出：最大的特征值 max
Begin
while (│diff│>ε) do
(1)for i=1 to n do
(1.1)sum=0
(1.2)for j= 1 to n do sum=sum+a[i,j]*x[j] end for

(1.3)b[i]= sum
end for (2)max=│b[1]│ (3)for i=2 to n do
if (│b[i]│>max) then max=│b[i]│ end if end for
(4)for i=1 to n do
x[i] =b[i]/max
end for
(5)diff=max-oldmax , oldmax=max
end while
End

1.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 维结果 向量中的最大值 localmax，然后通过归约操作的求最大值运算得到整个 n 维结果向量中的最 大值 max 并广播给所有处理器，各处理器利用 max 进行规格化操作。最后通过扩展收集操 作将各处理器中的 m 维结果向量按处理器编号连接起来并广播给所有处理器，以进行下一 次迭代计算，直至收敛。具体算法框架描述如下：
算法 21.2 乘幂法求解矩阵最大特征值的并行算法
输入：系数矩阵An×n ，初始向量v n×1 ，ε
输出：最大的特征值 max
Begin
对所有处理器 my_rank(my_rank=0,…, p-1)同时执行如下的算法:
while (│diff│>ε) do /* diff 为特征向量的各个分量新旧值之差的最大值*/ (1)for i=0 to m-1 do /*对所存的行计算特征向量的相应分量*/
(1.1)sum=0
(1.2)for j= 0 to n-1 do
sum=sum+a[i,j]*x[j]
end for
(1.3)b[i]= sum
end for
(2)localmax=│b[0]│ /*对所计算的特征向量的相应分量 求新旧值之差的最大值 localmax */
(3)for i=1 to m-1 do
if (│b[i]│>localmax) then localmax=│b[i]│ end if end for
(4)用 Allreduce 操作求出所有处理器中 localmax 值的最大值 max
并广播到所有处理器中
(5)for i=0to m-1 do /*对所计算的特征向量归一化 */

b[i] =b[i]/max
end for
(6)用 Allgather 操作将各个处理器中计算出的特征向量的分量的新值组合并广播到 所有处理器中
(7)diff=max-oldmax, oldmax=max
end 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阶实对称矩阵，则存在一个正交矩阵U，使得
UTAU=D

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


él1
ê
D = ê0
êM
ê
êë0


0 L
l2 L
M O
0 L


0 ù
ú
0 ú
M ú
ú
ln úû

这里λ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 = app cos2θ+aqq sin2θ+apq sin2θ bqq = app sin2θ+aqq cos2θ-apq sin2θ
bpq = ((aqq -app )sin2θ)/2+apq cos2θ bqp = bpq
bpj = apj cosθ+aqj sinθ bqj = -apj sinθ +aqj cosθ bip = aip cosθ+aiq sinθ biq = -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，则可令

qq
m = - a pq ,
n = 1 ( a

- a pp ) ,
w = sgn( n ) m

容易推出：


sin 2q = w ,
2

sin q =


2(1 +


w ,
1 - w 2 )


cosq =
m 2

1 - sin 2 q
+ n 2

利用sin2θ，sinθ和cosθ的值，即得矩阵B的各元素。矩阵A经过旋转变换，选定的非主

pq
对角元素apq 及aqp（一般是绝对值最大的）就被消去，且其主对角元素的平方和增加了 2a 2 ，

而非主对角元素的平方和减少了 2a2 pq ，矩阵元素总的平方和不变。通过反复选取主元素apq ， 并作旋转变换，就逐步将矩阵A变为对角矩阵。在对称矩阵中共有(n2-n)/2 个非主对角元素要 被消去, 而每消去一个非主对角元素需要对 2n个元素进行旋转变换, 对一个元素进行旋转 变换需要 2 次乘法和 1 次加法，若各取一次乘法运算时间或一次加法运算时间为一个单位时 间，则消去一个非主对角元素需要 3 个单位运算时间，所以下述算法 21.3 的一轮计算时间 复杂度为 (n2-n)/2*2n*3=3n2(n-1)=On3)。
算法 21.3 单处理器上雅可比法求对称矩阵特征值的算法
输入：矩阵An×n ，ε
输出：矩阵特征值 Eigenvalue
Begin
(1)while (max >ε) do (1.1) max=a[1,2] (1.2)for i=1 to n do
for j= i+1 to n do
if (│a[i,j]) │>max) then max =│a[i,j] │，p=i，q=j end if end for
end for
(1.3)Compute:
m=- a[p,q]，n=(a[q,q]- a[p,p])/2，w=sgn(n)*m/sqrt(m*m+nn)， 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)for j=1 to n do
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
end if end for
(1.5)for i=1 to n do

End
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]*co1
end if end for
(1.6)for i=1 to n do for j=1 to n do
a[i,j]=b[i,j]
end for end for
end while
(2)for i=1 to n do
Eigenvalue[i]= a[i, i]
end for

1.2.2 雅可比法求对称矩阵特征值的并行算法

串行雅可比算法逐次寻找非主对角元绝对值最大的元素的方法并不适合于并行计算。因 此，在并行处理中，我们每次并不寻找绝对值最大的非主对角元消去，而是按一定的顺序将 A 中的所有上三角元素全部消去一遍，这样的过程称为一轮。由于对称性，在一轮中，A 的 下三角元素也同时被消去一遍。经过若干轮，可使 A 的非主对角元的绝对值减少，收敛于 一个对角方阵。具体算法如下：
设处理器个数为p，对矩阵A按行划分为 2p块，每块含有连续的m/2 行向量，记 m = én / pù ， 这些行块依次记为A0 ,A1 , …,A2p-1 ，并将A2i 与A2i+1 存放与标号为i的处理器中。
每轮计算开始，各处理器首先对其局部存储器中所存的两个行块的所有行两两配对进行 旋转变换，消去相应的非对角线元素。然后按图 21.1 所示规律将数据块在不同处理器之间

传送，以消去其它非主对角元素。

图 1.1 p=4 时的雅可比算法求对称矩阵特征值的数据交换图

这里长方形框中两个方格内的整数被看成是所移动的行块的编号。在要构成新的行块配 对时，只要将每个处理器所对应的行块按箭头方向移至相邻的处理器即可，这样的传送可以 在行块之间实现完全配对。当编号为i和j的两个行块被送至同一处理器时，令编号为i的行块 中的每行元素和编号为j的行块中的每行元素配对，以消去相应的非主对角元素，这样在所 有的行块都两两配对之后，可以将所有的非主对角元素消去一遍，从而完成一轮计算。由图
1.1可以看出，在一轮计算中，处理器之间要 2p-1 次交换行块。为保证不同行块配对时可以 将原矩阵A的非主对角元素消去，引入变量b记录每个处理器中的行块数据在原矩阵A中的实 际行号。在行块交换时，变量b也跟着相应的行块在各处理器之间传送。
处理器之间的数据块交换存在如下规律：
0 号处理器前一个行块(简称前数据块，后一个行块简称后数据块)始终保持不变，将后 数据块发送给 1 号处理器，作为 1 号处理器的前数据块。同时接收 1 号处理器发送的后数据 块作为自己的后数据块。
p-1 号处理器首先将后数据块发送给编号为 p-2 的处理器，作为 p-2 号处理器的后数据 块。然后将前数据块移至后数据块的位置上，最后接收 p-2 号处理器发送的前数据块作为自 己的前数据块。
编号为 my_rank 的其余处理器将前数据块发送给编号为 my_rank+1 的处理器，作为 my_rank+1 号处理器的前数据块。将后数据块发送给编号为 my_rank-1 的处理器，作为 my_rank-1 号处理器的后数据块。
为了避免在通信过程中发生死锁，奇数号处理器和偶数号处理器的收发顺序被错开。
使偶数号处理器先发送后接收；而奇数号处理器先将数据保存在缓冲区中，然后接收偶数号 处理器发送的数据，最后再将缓冲区中的数据发送给偶数号处理器。数据块传送的具体过程 描述如下：
算法 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 号处理器发送来的后数据块作为自己的后数据块
end if
(2)if ((my-rank=p-1) and ( my-rank mod 2 ≠ 0)) then /*处理器 p-1 且其为奇数*/ (2.1)for i=m/2 to m-1 do /* 将后数据块保存在缓冲区 buffer 中*/
for j=0 to n-1 do
buffer[i-m/2,j]=a[i,j]
end for end for
(2.2)for i=m/2 to m-1 do

buf [i-m/2] =b[i]
end for
(2.3)for i=0 to m/2-1 do /*将前数据块移至后数据块的位置上*/
for j=0 to n-1 do
a[i+m/2,j]=a[i,j]
end for end for
(2.4)for i=0 to m/2-1 do
b[i+m/2] =b[i]
end for
(2.5)接收 p-2 号处理器发送的数据块作为自己的前数据块
(2.6)将 buffer 中的后数据块发送给编号为 p-2 的处理器
end if
(3)if ((my-rank=p-1) and ( my-rank mod 2=0)) then /*处理器 p-1 且其为偶数*/ (3.1)将后数据块发送给编号为 p-2 的处理器
(3.2)for i=0 to m/2-1 do /*将前数据块移至后数据块的位置上*/
for j=0 to n-1 do
a[i+m/2,j]=a[i,j]
end for end for
(3.3)for i=0 to m/2-1 do
b[i+m/2] =b[i]
end for
(3.4)接收 p-2 号处理器发送的数据块作为自己的前数据块
end if
(4)if ((my-rank ≠ p-1) and ( my-rank ≠ 0)) then /*其它的处理器*/ (4.1)if (my-rank mod 2=0) then /*偶数号处理器*/
(i)将前数据块发送给编号为 my_rank+1 的处理器
(ii)将后数据块发送给编号为 my_rank-1 的处理器
(ii)接收编号为 my_rank-1 的处理器发送的数据块作为自己的前 数据块
(iv)接收编号为 my_rank+1 的处理器发送的数据块作为自己的后 数据块
else /*奇数号处理器*/
(v)for i=0 to m-1 do /* 将前后数据块保存在缓冲区 buffer 中*/
for j=0 to n-1 do
buffer[i,j]=a[i,j]
end for end for
(vi)for i=0 to m-1 do
buf[i] =b[i]
end for
(vii)接收编号为 my_rank-1 的处理器发送的数据块作为自己的前 数据块

End


end if
(viii)接收编号为 my_rank+1 的处理器发送的数据块作为自己的 后数据块
(ix)将存于 buffer 中的前数据块发送给编号为 my_rank+1 的处理 器
(x)将存于 buffer 中的后数据块发送给编号为 my_rank-1 的处理器
end if

各处理器并行地对其局部存储器中的非主对角元素aij 进行消去，首先计算旋转参数并对 第i行和第j行两行元素进行旋转行变换。然后通过扩展收集操作将相应的旋转参数及第i列和 第j列的列号按处理器编号连接起来并广播给所有处理器。各处理器在收到这些旋转参数和 列号之后，按 0,1,…,p-1 的顺序依次读取旋转参数及列号并对其m行中的第i列和第j列元素进 行旋转列变换。
经过一轮计算的 2p-1 次数据交换之后，原矩阵 A 的所有非主对角元素都被消去一次。
此时，各处理器求其局部存储器中的非主对角元素的最大元 localmax，然后通过归约操作的 求最大值运算求得将整个 n 阶矩阵非主对角元素的最大元 max，并广播给所有处理器以决定 是否进行下一轮迭代。具体算法框架描述如下：
算法 21.5 雅可比法求对称矩阵特征值的并行算法
输入：矩阵An×n ，ε
输出：矩阵 A 的主对角元素即为 A 的特征值
Begin
对所有处理器 my_rank(my_rank=0,…, p-1)同时执行如下的算法:
(a)for i=0 to m-1 do
b[i] =myrank*m+i /* b 记录处理器中的行块数据在原矩阵 A 中的实际行号*/
end for
(b)while (│max│>ε) do /* max 为 A 中所有非对角元最大的绝对值*/ (1)for i=my_rank*m to (my_rank+1)*m-2 do
/*对本处理器内部所有行两两配对进行旋转变换*/
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 为它们在块内的相对行号*/
(1.2)if (a[r,j] ≠ 0) then /*对四个主元素的旋转变换*/ (i)Compute:
f=-a[r,j], g=( a[t,j]- a[r,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)，
bpp= a[r,i]*cos1*cos1+a[t,j]*sin1*sin1+a[r,j]*sin2， bqq= a[r,i]* sin1*sin1+a[t,j]* cos1*cos1-a[r,j]*sin2， bpq=0 , bqp=0
(ii)for v=0 to n-1 do /*对两个主行其余元素的旋转变换*/
if ((v ≠ i) and ( v ≠ j)) then
br[v] = a[r,v]*cos1 + a[t,v]*sin1
a[t,v] = -a[r,v]* sin1 + a[t,v]* cos1
end if

end for
(iii)for v=0 to n-1 do
if ((v ≠ i) and ( v ≠ j)) then
a[r,v]=br[v]
end if end for
(iv)for v=0 to m-1 do
/*对两个主列在本处理器内的其余元素的旋转变换*/
if (( v ≠ r) and ( v ≠ t)) then
bi[v] = a[v,i]*cos1 + a[v,j]*sin1
a[v,j]= - a[v,i]* sin1 + a[v,j]* cos1
end if end for
(v)for v=0 to m-1do
if ((v ≠ r) and ( v ≠ t)) then
a[v,i]= bi[v]
end if end for
(vi)Compute:
a[r,i]=bpp , a[t,j]=bqq , a[r,j]=bpq , a[t,i]=bqp，
/*用 temp1 保存本处理器主行的行号和旋转参数*/
temp1[0]=sin1, temp1[1]=cos1，
temp1[2]=(float)i ,temp1[3]= (float)j
else
(vii)Compute:
temp1[0]=0, temp1[1]= 0，
temp1[2]= 0 , temp1[3]= 0
end if
(1.3)将所有处理器 temp1 中的旋转参数及主行的行号 按处理器编号连接起来并广播给所有处理器,存于 temp2 中
(1.4)current=0
(1.5) for v=1 to p do
/*根据 temp2 中的其它处理器的旋转参数及主行的行号对相关的 列在本处理器的部分进行旋转变换*/
(i)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] (ii )if (s1、c1、 i1、 j1 中有一不为 0) then
if (my-rank ≠ current) then for z=0 to m-1 do
zi[z]=a[z,i1]*c1 + a[z,j1]*s1
a[ z,j1]=- a[z,i1]*s1 + a[z,j1]*c1
end for
for z=0 to m-1 do

a[zi1]= zi[z]
end for end if
end if
(iii)current=current+1
end for end for
end for
(2)for counter=1 to 2p-2 do
/*进行 2p-2 次处理器间的数据交换, 并对交换后处理器中所有行两两配对 进行旋转变换*/
(2.1)Data_exchange( ) /*处理器间的数据交换*/
(2.2)for i=0 to m/2-1 do
for j=m/2 to m-1 do
(i) if (a[i,b[j]] ≠ 0) then /*对四个主元素的旋转变换*/
①Compute:
f= -a[i,b[j]]，g=(a[j,b[j]]- a[i,b[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)，
bpp= a[i,b[i]]*cos1*cos1+ a[j,b[j]]*sin1*sin1+a[i,b[j]]*sin2， bqq= a[i,b[i]]* sin1*sin1+a[j,b[j]]* cos1*cos1-a[i,b[j]]*sin2， bpq=0, bqp=0
②for v=0 to n-1 do /*对两个主行其余元素的旋转变换*/
if ((v ≠ b[i]) and ( v ≠ b[j])) then
br[v] = a[i,v]*cos1 + a[j,v]*sin1
a[j,v] = -a[i,v]* sin1 + a[j,v]* cos1
end if end for
③for v=0 to n-1 do
if ((v ≠ b[i]) and ( v ≠ b[j])) then
a[i,v]=br[v]
end if end for
④for v=0 to m-1 do
/*对本处理器内两个主列的其余元素旋转变换*/
if ((v ≠ i) and ( v ≠ j)) then
bi[v] = a[v, b[i]]*cos1 + a[v, b[j]]*sin1
a[v, b[j]] = - a[v, b[i]]* sin1 + a[v, b[j]]* cos1
end if end for
⑤for v=0 to m-1 do
if ((v ≠ i) and ( v ≠ j)) then
a[v, b[i]]=bi[v]

end if end for
⑥Compute:
a[i, b[i]]=bpp , a[j, b[j]]=bqq ,
a[i, b[j]]=bpq , a[j, b[i]]=bqp
/*用 temp1 保存本处理器主行的行号和旋转参数*/
temp1[0]=sin1 , temp1[1]=cos1,
temp1[2]=(float)b[i] , temp1[3]= (float)b[j]

end for
else
⑦Compute: temp1[0]=0,temp1[1]= 0, temp1[2]= 0,temp1[3]= 0
end if
(ii)将所有处理器 temp1 中的旋转参数及主行的行号按处理 器编号连接起来并广播给所有处理器,存于 temp2 中
(iii)current=0 (iv)for v=1 to p do
/*根据 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) then if (my-rank ≠ current) then
for z=0 to m-1 do
zi[z]=a[z,i1]*c1 + a[z,j1]*s1
a[ z,j1]=- a[z,i1]s1 + a[z,j1]*c1
end for
for z=0 to m-1 do
a[z,i1]= zi[z]
end for end if
end if
③current=current+1
end for end for
end for

(3)Data-exchange( )
/*进行一轮中的最后一次处理器间的数据交换,使数据回到原来的位置*/ (4)localmax=max /*localmax 为本处理器中非对角元最大的绝对值/
(5)for i=0 to m-1 do for j=0 to n-1 do
if ((m*my-rank+i) ≠ j) then

End
if (│a[i,j]│> localmax) then localmax=│a[i,j]│ end if end if
end for end for
(6)通过 Allreduce 操作求出所有处理器中 localmax 的最大值为新的 max 值
end while

在上述算法中, 每个处理器在一轮迭代中要处理 2 p个行块对, 由于每一个行块对含有 m/2 行, 因而对每一个行块对的处理要有(m/2)2=m2/4 个行的配对, 即消去m2/4 个非主对角 元素. 每消去一个非主对角元素要对同行n 个元素和同列n 个元素进行旋转变换. 由于按行 划分数据块, 对同行n 个元素进行旋转变换的过程是在本处理器中进行的. 而对同列n 个元 素进行旋转变换的过程则分布在所有处理器中进行. 但由于所有处理器同时在为其它处理 器的元素消去过程进行列的旋转变换，对每个处理器而言，对列元素进行旋转变换的处理总 量仍然是n 个元素。 因此，一个处理器每消去一个非主对角元素共要对 2n 个元素进行旋 转变换。而对一个元素进行旋转变换需要 2 次乘法和 1 次加法，若取一次乘法运算时间或一 次加法运算时间为一个单位时间，则其需要 3 个单位运算时间，即一轮迭代的计算时间为T1
＝3×2 p ×2n ×m2/4＝3n3/p；在每轮迭代中，各个处理器之间以点对点的通信方式相互错
开交换数据 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 + 2]t s (
p - 1) + [2n(n - 1) / p + 1]t w ( p - 1)

因此雅可比算法求对矩阵特征值的一轮并行计算时间为Tp = T1 + T2 。我们的大量实验结果

说明，对于n阶方阵，用上述算法进行并行计算，一般需要不超过O(logn)轮就可以收敛。
MPI 源程序请参见章末附录。

1.3 求对称矩阵特征值的单侧旋转法

1.3.1 单侧旋转法的算法描述
求解对称方阵特征值及特征向量的雅可比算法在每次消去计算前都要对非主对角元素 选择最大元，导致非实际计算开销很大。在消去计算时，必须对方阵同时进行行、列旋转变 换，这称之为双侧旋转(Two-side Rotation)变换。在双侧旋转变换中，方阵的行、列方向都 有数据相关关系，使得整个计算中的相关关系复杂，特别是在对高阶矩阵进行特征值求解时， 增加了处理器间的通信开销。针对传统雅可比算法的缺点，可以将双侧旋转改为单侧旋转， 得出一种求对称矩阵特征值的快速算法。这一算法对矩阵仅实施列变换，使得数据相关关系 仅在同列之间，因此方便数据划分，适合并行计算，称为单侧旋转法（One-side Rotation）。 若A 为一对称方阵， λ 是A 的特征值，x 是A 的特征向量，则有：Ax=λx 。又A=AT ，所以 ATAx=Aλx=λλx，这说明了 λ2是ATA的特征值，x是ATA的特征向量 ，即ATA的特征值是A的特 征值的平方，并且它们的特征向量相同。
我们使用 18.7.1 节中所介绍的Givens旋转变换对A进行一系列的列变换，得到方阵Q使 其各列两两正交，即AV＝Q，这里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 =δ ,

i i i
i i i

i
||qi ||2 =
qT qi =
d i = li
。因此在将A进行列变换得到各列两两正交的方阵Q后，其各列的

谱范数||qi ||2 即为A的特征值的绝对值。记V的第i列为vi ， AV＝Q，所以Avi = qi 。又因为vi 为 A的特征向量，所以Avi =λi vi ，即推得λi vi = qi 。因此λi 的符号可由向量 qi 及vi 的对应分量是 否同号判别，实际上在具体算法中我们只要判断它们的第一个分量是否同号即可。若相同， 则λi 取正值，否则取负值。
求对称矩阵特征值的单侧旋转法的串行算法如下： 算法 21.6 求对称矩阵特征值的单侧旋转法
输入：对称矩阵 A，精度 e
输出：矩阵 A 的特征值存于向量 B 中
Begin
(1)while (p>ε)
p=0
for i=1 to n do
for j=i+1 to n do
(1.1 )sum[0]=0, sum[1]=0, sum[2]=0 (1.2 )for k=1 to n do
sum[0]= sum[0]+a[k,i]*a[k, j] sum[1]= sum[1]+ a[k,i]* a[k,i] sum[2]= sum[2]+ a[k, j]* a[k, j]
end for
(1.3 )if (│sum[0]│>ε) then
(i) aa=2*sum[0] bb=sum[1]-sum[2]