9矩阵特征值计算 在实际的工程计算中,经常会遇到求n阶方阵A的特征值( Eigenvalue)与特征向量 ( Eigenvector)的问题。对于一个方阵A,如果数值λ使方程组 即(4-)x=0有非零解向量 Solution vector)x,则称λ为方阵A的特征值,而非零向量x为 特征值所对应的特征向量,其中Ln为n阶单位矩阵 由于根据定义直接求矩阵特征值的过程比较复杂,因此在实际计算中,往往采取一些 数值方法。本章主要介绍求一般方阵绝对值最大的特征值的乘幂( Power)法、求对称方阵特 征值的雅可比法和单侧旋转(One- side rotation)法以及求一般矩阵全部特征值的QR方法及 一些相关的并行算法 1.1求解矩阵最大特征值的乘幂法 111乘幂法及其串行算法 在许多实际问题中,只需要计算绝对值最大的特征值,而并不需要求矩阵的全部特征值 乘幂法是一种求矩阵绝对值最大的特征值的方法。记实方阵A的n个特征值为A (1,2,…,m),且满足: 2|2|≥|a32…≥|n 特征值λ对应的特征向量为x。乘幂法的做法是:①取n维非零向量w作为初始向量:② 对于k=1 ,做如下迭代: uk=Avk- V=Wk/uk 直至脚b-ke)de ()for i=l to n do (1.2)for=I to n do sm=sm+a[i引*xU enal io
1. 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, …,做如下迭代： uk =Avk-1 vk= uk /║uk║∞ 直至 − uk+1 uk ε 为止，这时 vk+1 就是 A 的绝对值最大的特征值 λ1 所对应的特 征向量 x1。若 vk-1 与 vk的各个分量同号且成比例，则 λ1=║uk║∞；若 vk-1 与 vk的各个分量异号 且成比例，则 λ1= -║uk║∞。若各取一次乘法和加法运算时间、一次除法运算时间、一次比较 运算时间为一个单位时间，则因为一轮计算要做一次矩阵向量相乘、一次求最大元操作和一 次规格化操作，所以下述乘幂法串行算法 21.1 的一轮计算时间为 n 2+2n=O(n 2 )。 算法 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}=sm (2)max b[ (3)for i=2 to n d if( b(a >max)then max= b(i end if end fo (4)for i=I to n do x[=b/ end for (5)diff=max-oldmax, oldmaxmax end while 112乘幕法并行算法 乘幂法求矩阵特征值由反复进行矩阵向量相乘来实现,因而可以采用矩阵向量相乘的数 据划分方法。设处理器个数为p,对矩阵A按行划分为p块,每块含有连续的m行向量 这里m=「m/pl,编号为i的处理器含有A的第m至第(+1)m-1行数据,(=01,…P-1),初 始向量γ被广播给所有处理器 各处理器并行地对存于局部存储器中a的行块和向量ν做乘积操作,并求其m维结果 向量中的最大值 localmax,然后通过归约操作的求最大值运算得到整个n维结果向量中的最 大值max并广播给所有处理器,各处理器利用max进行规格化操作。最后通过扩展收集操 作将各处理器中的m维结果向量按处理器编号连接起来并广播给所有处理器,以进行下 次迭代计算,直至收敛。具体算法框架描述如下: 算法212乘幂法求解矩阵最大特征值的并行算法 输入:系数矩阵Amxn,初始向量νn×1,ε 输出:最大的特征值max Begin 对所有处理器 my rank( my rank=0,…,p-1)同时执行如下的算法 hil(dile)do/dif.特征向量的各个分量新旧值之差的最大值* (1)for=0tom-1do/*对所存的行计算特征向量的相应分量* (1.1)am=0 (1.2)forj=0 to n-l do sumsumtai, xlI end fo (1.3)6[}= end for (2) calmar=|b0]|/对所计算的特征向量的相应分量 求新旧值之差的最大值 locaImax+ ()for i=l to m-I do r( b[>localmax)then localmax b[a end if (4)用 Allreduce操作求出所有处理器中 locaImax值的最大值max 并广播到所有处理器中 (5)fori=0tom-1do/*对所计算的特征向量归一化*
(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 /*对所计算的特征向量归一化 */
b0=b[amax end for 6)用A! gather操作将各个处理器中计算出的特征向量的分量的新值组合并广播到 所有处理器中 (7)diff=max-oldmax, oldmaxmax end while 若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时 间,一轮迭代的计算时间为m+2m;一轮迭代中,各处理器做一次归约操作,通信量为1, 次扩展收集操作,通信量为m,则通信时间为4,(√p-1)+(m+1)n(p-1)。因此乘幂法的 轮并行计算时间为Tn=m+2m+4(√P-1)+(m+1xn(P-1) MPI源程序请参见所附光盘。 12求对称矩阵特征值的雅可比法 121雅可比法求对称矩阵特征值的串行算法 若矩阵A=[a是n阶实对称矩阵,则存在一个正交矩阵U,使得 UAU=D 其中D是一个对角矩阵,即 0 这里λ(=1,2…,m)是A的特征值,U的第i列是与λ对应的特征向量。雅可比算法主要是通 过正交相似变换将一个实对称矩阵对角化,从而求出该矩阵的全部特征值和对应的特征向 量。因此可以用一系列的初等正交变换逐步消去A的非对角线元素,从而使矩阵A对角化。 设初等正交矩阵为R(PqO,其(pP)(q,q)位置的数据是cosO,(p,q)(qp)位置的数据分别是 sinO和 sinA(p≠q),其它位置的数据和一个同阶数的单位矩阵相同。显然可以得到 不妨记B=R(PqO)AR(PqO,如果将右端展开,则可知矩阵B的元素与矩阵A的元素 之间有如下关系 bpp=appcos20+aggsin20+apgsin20 bpg=((aqg-app)sin20)/2+apg cos bpi=apjcoso+aa sine bip=aip cose+ aigsine bi=ajr 其中1≤,j≤n且≠Pq。因为A为对称矩阵,R为正交矩阵,所以B亦为对称矩阵 若要求矩阵B的元素b=0,则只需令(aq-ap)sin2)2+apqc0s20=0,即:
b[i] =b[i]/max end for (6)用 Allgather 操作将各个处理器中计算出的特征向量的分量的新值组合并广播到 所有处理器中 (7)diff=max-oldmax, oldmax=max end while End 若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时 间，一轮迭代的计算时间为 mn+2m；一轮迭代中，各处理器做一次归约操作，通信量为 1， 一次扩展收集操作，通信量为 m，则通信时间为 4t ( p −1) + (m +1)tw( p −1) s 。因此乘幂法的 一轮并行计算时间为 T = mn + 2m + 4t ( p −1) + (m +1)tw(p −1) p s 。 MPI 源程序请参见所附光盘。 1.2 求对称矩阵特征值的雅可比法 1.2.1 雅可比法求对称矩阵特征值的串行算法 若矩阵 A=[aij]是 n 阶实对称矩阵，则存在一个正交矩阵 U，使得 U TAU=D 其中 D 是一个对角矩阵，即 = n D 0 0 0 0 0 0 2 1 这里 λ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θ+aqq cos2θ-apqsin2θ bpq = ((aqq -app)sin2θ)/2+apq cos2θ 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+apq cos2θ=0，即：
1g20=(a-am)/2 在实际应用时,考虑到并不需要解出O,而只需要求出sin26,sin0和co就可以了。 如果限制θ值在-x/2e)de (1.1)max=a[12] (1.2)for i=l to n de if(a[/)>max)then max=a[], P=i, qj end if end for m-alp, q ], n(alg. q-alp pl)2, w=sgn(n)*m/sqrt(m*m+n*n si2=w,sil=w/sqrt(2(1+ sqrt(l-l*w))), col=sqrt( blp,P=alp.]*col *col+al, q*si l*sil+ alp. q]*si2 bl, q=alp,P]*sil*sil+ alg, q*col*col-alp, q*si2 b{q,p]=0,b{P,q}=0 if(≠p)and(≠q))then (ibp=alp col+ala sil (1. 5)for i=l to n do
( ) 2 2 qq pp pq a a a tg − − = 在实际应用时，考虑到并不需要解出 θ，而只需要求出 sin2θ，sinθ 和 cosθ 就可以了。 如果限制 θ 值在-π/2 ε) 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+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)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
if(i+p)and(i≠q)then (ib, P=aip]col+ai, q* sil (ii)b[i, g=-alip*sil+ ali, q]"col end if end for (1. 6for i=/ to n do forj=l to n de q[/=b end for end while (2)for i=l to n do Eigenvalue=a[i, i] end for End 122雅可比法求对称矩阵特征值的并行算法 串行雅可比算法逐次寻找非主对角元绝对值最大的元素的方法并不适合于并行计算。因 此,在并行处理中,我们每次并不寻找绝对值最大的非主对角元消去,而是按一定的顺序将 A中的所有上三角元素全部消去一遍,这样的过程称为一轮。由于对称性,在一轮中,A的 下三角元素也同时被消去一遍。经过若干轮,可使A的非主对角元的绝对值减少,收敛于 个对角方阵。具体算法如下 设处理器个数为p,对矩阵A按行划分为2p块,每块含有连续的m/2行向量,记 m=「n/p1,这些行块依次记为A4,…Ayp,并将A2与A2+存放与标号为i的处理器中 每轮计算开始,各处理器首先对其局部存储器中所存的两个行块的所有行两两配对进行 旋转变换,消去相应的非对角线元素。然后按图21.1所示规律将数据块在不同处理器之间 传送,以消去其它非主对角元素。 开始:(0,1)(2,3)(4,5)(6,7) 七步: 图11p=4时的雅可比算法求对称矩阵特征值的数据交换图
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 End 1.2.2 雅可比法求对称矩阵特征值的并行算法 串行雅可比算法逐次寻找非主对角元绝对值最大的元素的方法并不适合于并行计算。因 此，在并行处理中，我们每次并不寻找绝对值最大的非主对角元消去，而是按一定的顺序将 A 中的所有上三角元素全部消去一遍，这样的过程称为一轮。由于对称性，在一轮中，A 的 下三角元素也同时被消去一遍。经过若干轮，可使 A 的非主对角元的绝对值减少，收敛于 一个对角方阵。具体算法如下： 设处理器个数为 p，对矩阵 A 按行划分为 2p 块，每块含有连续的 m/2 行向量，记 m = n / p ，这些行块依次记为 A0,A1, …,A2p-1，并将 A2i 与 A2i+1 存放与标号为 i 的处理器中。 每轮计算开始，各处理器首先对其局部存储器中所存的两个行块的所有行两两配对进行 旋转变换，消去相应的非对角线元素。然后按图 21.1 所示规律将数据块在不同处理器之间 传送，以消去其它非主对角元素。 开 始 ：（0，1）（2，3）（4，5）（6，7） 第一步 ：（0，3）（1，5）（2，7）（4，6） 第二步 ：（0，5）（3，7）（1，6）（2，4） 第三步 ：（0，7）（5，6）（3，4）（1，2） 第四步 ：（0，6）（7，4）（5，2）（3，1） 第五步 ：（0，4）（6，2）（7，1）（5，3） 第六步 ：（0，2）（4，1）（6，3）（7，5） 第七步 ：（0，1）（2，3）（4，5）（6，7） 0 1 2 3 4 5 6 7 图 1.1 p=4 时的雅可比算法求对称矩阵特征值的数据交换图
这里长方形框中两个方格内的整数被看成是所移动的行块的编号。在要构成新的行块配 对时,只要将每个处理器所对应的行块按箭头方向移至相邻的处理器即可,这样的传送可以 在行块之间实现完全配对。当编号为i和j的两个行块被送至同一处理器时,令编号为i的 行块中的每行元素和编号为j的行块中的每行元素配对,以消去相应的非主对角元素,这样 在所有的行块都两两配对之后,可以将所有的非主对角元素消去一遍,从而完成一轮计算 由图1.1可以看出,在一轮计算中,处理器之间要2p-1次交换行块。为保证不同行块配对 时可以将原矩阵A的非主对角元素消去,引入变量b记录每个处理器中的行块数据在原矩 阵A中的实际行号。在行块交换时,变量b也跟着相应的行块在各处理器之间传送。 处理器之间的数据块交换存在如下规律 0号处理器前一个行块(简称前数据块,后一个行块简称后数据块)始终保持不变,将后 数据块发送给1号处理器,作为1号处理器的前数据块。同时接收1号处理器发送的后数据 块作为自己的后数据块 1号处理器首先将后数据块发送给编号为p-2的处理器,作为p-2号处理器的后数据 块。然后将前数据块移至后数据块的位置上,最后接收p-2号处理器发送的前数据块作为自 己的前数据块。 编号为 my rank的其余处理器将前数据块发送给编号为 my rank+1的处理器,作为 my rank+Ⅰ号处理器的前数据块。将后数据块发送给编号为 my rank-1的处理器,作为 my rank-1号处理器的后数据块。 为了避免在通信过程中发生死锁,奇数号处理器和偶数号处理器的收发顺序被错开。 使偶数号处理器先发送后接收:而奇数号处理器先将数据保存在缓冲区中,然后接收偶数号 处理器发送的数据,最后再将缓冲区中的数据发送给偶数号处理器。数据块传送的具体过程 描述如下: 算法214雅可比法求对称矩阵特征值的并行过程中处理器之间的数据块交换算法 输入:矩阵A的行数据块和向量b的数据块分布于各个处理器中 输出:在处理器阵列中传送后的矩阵A的行数据块和向量b的数据块 对所有处理器 my rank( my rank=0…,p-1)同时执行如下的算法 矩阵A和向量b为要传送的数据块* (1)if( my-rank=0)then/*0号处理器* (1.1)将后数据块发送给1号处理器 (1.2)接收1号处理器发送来的后数据块作为自己的后数据块 (2)if( my-rank=p-1)and(my- rank mod2≠0)then/*处理器p-1且其为奇数* (2.1)for=m/2tom-ldo/*将后数据块保存在缓冲区buer中* forj=0 to n-l do buffer(i-m/2=a[ij d fo 2.2 )for i=m/2 to m-I do
这里长方形框中两个方格内的整数被看成是所移动的行块的编号。在要构成新的行块配 对时，只要将每个处理器所对应的行块按箭头方向移至相邻的处理器即可，这样的传送可以 在行块之间实现完全配对。当编号为 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]=bd (2.3)for=0tom2-ldo/*将前数据块移至后数据块的位置上*/ for产=0ton-ldo a计+m2=a[i end for (2. 4)for i=0 to m/2-1 do b计+m/2]=b d fo (2.5)接收p-2号处理器发送的数据块作为自己的前数据块 (26)将bfer中的后数据块发送给编号为p2的处理器 3)if( my-rank=p1)and( my-rank mod2=0)then/*处理器p1且其为偶数* (3.1)将后数据块发送给编号为p2的处理器 (3,2)fori=0tom/2-1do/*将前数据块移至后数据块的位置上* orj=0 to n-l do a计+m2=a nd fo (3. 3)for i=0 to m/2-1 do b计+m/2]=b end ior (34)接收p-2号处理器发送的数据块作为自己的前数据块 (4)if(( my-rank≠p-l)and( my-rank≠0)then/其它的处理器* (4ji( my-rank mod2=0)then/偶数号处理器*/ ()将前数据块发送给编号为 my rank+1的处理器 (i)将后数据块发送给编号为 my rank-1的处理器 (in)接收编号为 my rank-1的处理器发送的数据块作为自己的前 数据块 (iv)接收编号为 my rank+1的处理器发送的数据块作为自己的后 数据块 else/*奇数号处理器* (v)for=0tom-1do/*将前后数据块保存在缓冲区buer中* for户=0ton-ldo buffer=a[i, end for end for (vifor i=0 to m-l do buf[]=b[i end for (ⅶi)接收编号为 my rank-1的处理器发送的数据块作为自己的前 数据块
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 的处理器发送的数据块作为自己的前 数据块
(vi)接收编号为 my rank+1的处理器发送的数据块作为自己的 后数据块 (ⅸx)将存于buer中的前数据块发送给编号为 my rank+1的处理 (x)将存于buer中的后数据块发送给编号为 my rank-1的处理器 end if end if End 各处理器并行地对其局部存储器中的非主对角元素a进行消去,首先计算旋转参数并 对第i行和第j行两行元素进行旋转行变换。然后通过扩展收集操作将相应的旋转参数及第 i列和第j列的列号按处理器编号连接起来并广播给所有处理器。各处理器在收到这些旋转 参数和列号之后,按0,1,…P-1的顺序依次读取旋转参数及列号并对其m行中的第i列和第 j列元素进行旋转列变换 经过一轮计算的2p-1次数据交换之后,原矩阵A的所有非主对角元素都被消去一次 此时,各处理器求其局部存储器中的非主对角元素的最大元 locaImax,然后通过归约操作的 求最大值运算求得将整个n阶矩阵非主对角元素的最大元muax,并广播给所有处理器以决定 是否进行下一轮迭代。具体算法框架描述如下: 算法215雅可比法求对称矩阵特征值的并行算法 输入:矩阵Anxn,ε 输出:矩阵A的主对角元素即为A的特征值 对所有处理器 my rank( my rank=0,…,p-1)同时执行如下的算法 (afor i=0 to m-l do b[]= myrank*m+i/b记录处理器中的行块数据在原矩阵A中的实际行号* end for (b)whie(|max|>e)do/*max为A中所有非对角元最大的绝对值 (I)for i=my rank*m to(my rank+l)*m-2 do /*对本处理器内部所有行两两配对进行旋转变换* for j=i+l to(my rank+l)*m-1 do (1.1)= i mod m,闩modm陣*,j为进行旋转变换行(称为主行)的 实际行号,r,t为它们在块内的相对行号 (1.2ff(ar≠0)then/*对四个主元素的旋转变换 fa[rl,g=(a[小-a[,)/2,h=gm(g)°sqgr(+g*g), sin2=h, sinl=h/sgri(2 (1+sgrt(l-hh)), cosl=sqrt(l-sin1sinl) p=, i*cosl*cosl+at,*sinl"sinl+a[; /sin2 bear i" sinI sinl+at," cosl"cosl-arjsin2 (i)fory=0ton-1do体*对两个主行其余元素的旋转变换* if((v≠)and(v≠)then 付y]=a[r;v]*cosl+qy a[t, v]=-a[r, v]* sinl alt, vI
(viii)接收编号为 my_rank+1 的处理器发送的数据块作为自己的 后数据块 (ix)将存于 buffer 中的前数据块发送给编号为 my_rank+1 的处理 器 (x)将存于buffer中的后数据块发送给编号为my_rank-1的处理器 end if end if End 各处理器并行地对其局部存储器中的非主对角元素 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
(infor v=0 to n-l do if((v≠)and(v≠)then end if end for (iv)for I=0 to m-l do /*对两个主列在本处理器内的其余元素的旋转变换* if(v≠r)and(v≠)then biv]=av, i* av小]*sinl av=-alv i*sinl +alv* cosl if(v≠p)and(v≠D))then a[, i]= bilv] a[r, i]=bpp, alt=bqg, a[r=bpq, a[t, i]=bqp /*用 templ保存本处理器主行的行号和旋转参数* templ[0]=sinl, templ[l=cosl templ[2=(float)i, templ]=(float y (vii )Compute templ[0]=0, templ[1=0 templ[2]=0, templ3=0 end if (1.3)将所有处理器empl中的旋转参数及主行的行号 按处理器编号连接起来并广播给所有处理器,存于temp2中 (1.5)for I=l to p do *根据temp2中的其它处理器的旋转参数及主行的行号对相关的 列在本处理器的部分进行旋转变换* (i)Compute sl=temp2[(-1)*4+0],cl=temp2(y-1)*4+1], il=(int)temp2((1-1)*4+2],j1=(int )temp2[(1-1)*4+3] (i)if(sl、cl、il、jl中有一不为o)then if (my-rank current)then zi[-=a[=,1]'cl+a{=/1]*sl q[=/l]=a[=,]*sl+a[=/1]*cl end for for=0 to m- do
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=,il=-i{- end for nd if end if (iii)currentcurrent+I end f d fo (2)for counter=l to 2p-2 do *进行2p2次处理器间的数据交换,并对交换后处理器中所有行两两配对 进行旋转变换* (2.1) Data exchange()/处理器间的数据交换* (2.2)for i=0 to m/2-1 do for j=m/2 to m-I do (i)if(a[i,b]≠0)then/*对四个主元素的旋转变换* ①Com f-a[i, bll, g=(al, bll-ali, b[aD)/2 sgn(g)*sgp(升+g*g), sin2=h, sinl=h/sqrt(2*(1+sgr(1-h*h))) bpp=a[i, b[i*cosl*cosl+ a[, bll*sinl*sinl+a[i, b00*sin2 bgg=a[i, b[il* sinl*sinl+al, b l]*cosl *cosl-a[i, b0)*sin2 bpo, bgp=0 ②fory=0ton-1do/*对两个主行其余元素的旋转变换 if((v≠band(v≠bU])then brlv]=a[i, v*cosl +al, v]*sinl al, v]=-a[i, v]* sin1 +al, v]*cosl end if end for ③for= if(v≠b[)and(v≠b))then ali, vI=brlvl end if end for ④for1=0tom-1 /*对本处理器内两个主列的其余元素旋转变换 if(v≠)and(v≠))then biv]=av, b[i*cosl +av, blJ*sin a{v,b]]=-a[v,b[小*sin1+av,b]°cosl end for ⑤for=0tom-ldo if((v≠D)and(V≠)then Iv, ba=bilv
a[z,i1]= 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]