end if ⑥ a[i, b[0=bpp, al, blJ=bqq a[i, blll=bpg, al, b[al=bqp /*用empl保存本处理器主行的行号和旋转参数* mpl[o=sinl templ[2]=(float)b[i, templ3=(float)ol1 else templ[0=0, templ[l=0 end if (i)将所有处理器 templ中的旋转参数及主行的行号按处理 器编号连接起来并广播给所有处理器,存于lemp2中 (currer (iv )for I=l to p do /*根据temp2中的其它处理器的旋转参数及主行的行号 对相关的列在本处理器的部分进行旋转变换* sl=temp2[(y-1)*4+0],cl=temp2[(-1)*4+1] il-( int)temp2[(v-1)*4+2],=int)emp2[(v-1)*4+3] ②i(sl、cl、ilj1中有一不为0)then if( my-rank≠ current) then for=0 to m-l de -=a[=,i*cl+a[二/1]*sl a[=j1]=-a=, il]sl+ a= cl r=0 to m-1 do r,il==iEl d if nd fo endl ior (3)Data-exchange 进行一轮中的最后一次处理器间的数据交换,使数据回到原来的位置* (4) Dlocalmax=mx/* localmax为本处理器中非对角元最大的绝对值* ()for i=0 to m-l de for户=0ton-ld if(m*my-rank+)≠then
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] 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 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
if(ai/> localmax)then localmar= a[i,/ end if end if (6)通过 Allreduce操作求出所有处理器中 localmax的最大值为新的max值 End 在上述算法中,每个处理器在一轮迭代中要处理2p个行块对,由于每一个行块对含有 m/2行,因而对每一个行块对的处理要有(m/2)2=m24个行的配对,即消去m214个非主对角 元素.每消去一个非主对角元素要对同行n个元素和同列n个元素进行旋转变换由于按 行划分数据块,对同行n个元素进行旋转变换的过程是在本处理器中进行的.而对同列n 个元素进行旋转变换的过程则分布在所有处理器中进行但由于所有处理器同时在为其它 处理器的元素消去过程进行列的旋转变换,对每个处理器而言,对列元素进行旋转变换的处 理总量仍然是n个元素。因此,一个处理器每消去一个非主对角元素共要对2n个元素进 行旋转变换。而对一个元素进行旋转变换需要2次乘法和1次加法,若取一次乘法运算时间 或一次加法运算时间为一个单位时间,则其需要3个单位运算时间,即一轮迭代的计算时间 为71=3×2p×2nxm214=3n3lp;在每轮迭代中,各个处理器之间以点对点的通信方式相 互错开交换数据2p-1次,通信量为m+m,扩展收集操作n(n1)(2p)次,通信量为4,另外 有归约操作1次,通信量为1,从而不难得出上述算法求解过程中的总通信时间为: 7=[4,+m(n+1)n]4p-2)+[m(n-1)/p+2p(√p-1)+[2m(n-1)/p+1n(P-1) 因此雅可比算法求对矩阵特征值的一轮并行计算时间为T,=T+T2。我们的大量实验结果 说明,对于n阶方阵,用上述算法进行并行计算,一般需要不超过O(logn)轮就可以收敛。 MPⅠI源程序请参见章末附录 13求对称矩阵特征值的单侧旋转法 131单侧旋转法的算法描述 求解对称方阵特征值及特征向量的雅可比算法在每次消去计算前都要对非主对角元素 选择最大元,导致非实际计算开销很大。在消去计算时,必须对方阵同时进行行、列旋转变 换,这称之为双侧旋转(Two- side Rotation)变换。在双侧旋转变换中,方阵的行、列方向都 有数据相关关系,使得整个计算中的相关关系复杂,特别是在对高阶矩阵进行特征值求解时, 增加了处理器间的通信开销。针对传统雅可比算法的缺点,可以将双侧旋转改为单侧旋转, 得出一种求对称矩阵特征值的快速算法。这一算法对矩阵仅实施列变换,使得数据相关关系 仅在同列之间,因此方便数据划分,适合并行计算,称为单侧旋转法(One- side rotation) 若A为一对称方阵,λ是A的特征值,x是A的特征向量,则有:Ax=Ax。又A=A,所以 AAx=A入x=x,这说明了22是AA的特征值,x是AA的特征向量,即AA4的特征值是A 的特征值的平方,并且它们的特征向量相同 我们使用1871节中所介绍的 Givens旋转变换对A进行一系列的列变换,得到方阵Q 使其各列两两正交,即AV=Q,这里V为表示正交化过程的变换方阵。因QQ=(A)= AA=diag(b1,62,…,8n)为n阶对角方阵,可见这里61,b2,…,bn是矩阵AA的特征
if (│a[i,j]│> localmax) then localmax=│a[i,j]│ end if end if end for end for (6)通过 Allreduce 操作求出所有处理器中 localmax 的最大值为新的 max 值 end while End 在上述算法中, 每个处理器在一轮迭代中要处理 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＝3n 3 /p；在每轮迭代中，各个处理器之间以点对点的通信方式相 互错开交换数据 2p-1 次，通信量为 mn+m，扩展收集操作 n(n-1)/(2p)次，通信量为 4，另外 有归约操作 1 次，通信量为 1，从而不难得出上述算法求解过程中的总通信时间为： [4 ( 1) ](4 2) [ ( 1)/ 2] ( 1) [2 ( 1)/ 1] ( 1) T2 = t s + m n + tw p − + n n − p + t s p − + n n − p + tw p − 因此雅可比算法求对矩阵特征值的一轮并行计算时间为 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=A T，所以 A TAx=Aλx=λλx，这说明了 λ 2是 A TA 的特征值，x 是 A TA 的特征向量 ，即 A TA 的特征值是 A 的特征值的平方，并且它们的特征向量相同。 我们使用 18.7.1 节中所介绍的 Givens 旋转变换对 A 进行一系列的列变换，得到方阵 Q 使其各列两两正交，即 AV＝Q，这里 V 为表示正交化过程的变换方阵。 因 Q TQ=(AV) TAV=V TA TAV = diag(δ1,δ2, …,δn) 为 n 阶对角方阵，可见这里δ1,δ2, …,δn 是矩阵 A TA 的特征