求解大规模矩阵特征问题的并行算法研究

iter=iter+1;
6 if t < alpha*z and iter < 2, then go to 3
7 hj+1, j=sqrt(t);
8
; v
r j+1
=w
r j+1
/h
j+1,j
end
在算法 1 中，H(1:j, j)表示(h1, j, h2, j,…, hj, j)T；func(A, vrj )
数，若所有残量范数均小于允许误差，则停止迭代。否则转
第(5)步。
(5)计算 ur = λ1u1 + λ2u2 +" + λk / Puk / P ，, r=1,2,…,P，使用 全局归约求和函数得到 u=u1+u2+…+up，然后将 u 全局广播
σ min (H m − λi I) 的右奇异向量，那么有
ui = Vm zi
(
1
)
|| ( A − λi I )ui ||2 = σ min (H m − λi I)
(2)
其中， I 的前 m 行是 m 阶单位矩阵，最后一行的元素的值均
等于零。
该定理说明了求解精化向量只需只须计算(m+1)×m 阶矩 阵的奇异值分解。，并且用 σ min (H m − λi I) 判断 Ritz 值和精化
( λ , x)即可作为目标特征对的好的近似，从而停止迭代。否
则由 Ritz 向量构造 v1 进行重启动。
2.2 精化向量
精化重启动 Arnoldi 方法对于 Ritz 值 λi 利用满足如下条
件的精化向量 ui 作为目标特征向量的近似[2]：。
||
(A
−
λi I )ui
||2 =
min
u∈K (v1,A),||u||2 =1
素 hm+1,m 的值不等于零；em 是第 m 个坐标向量。
假设 Hm y = λ y ，则向量 x=Vmy 满足：
|| r ||2 =|| Ax − λ x ||2 =| hm+1,mem* y |
通常称 x 为 Ritz 向量， λ 为 Ritz 值。重启动 Arnoldi 算
法用 Ritz 对作为目标特征对的近似。如果||r||2 足够小，那么
基金项目：国家自然科学基金资助项目 “当代并行机的并行算法应 用基础研究”(60533020)；中国家科学院“十一五”信息化专项基金 资助项目“超级计算环境建设与应用”(INFO-115-B01) 作者简介：赵 韬(1979－)，男，博士研究生，主研方向：大规模非 对称矩阵特征问题的并行计算；迟学斌、陆忠华，研究员、博士生 导师；赵永华，副研究员 收稿日期：2009-10-17 E-mail：zhaotao@
器进程得到 Hm 和 Vm+1 相应的行块。 (3)若 hm+1,m=0，则计算 Hm 的全部特征值 λi 和 Ritz 向量
xi，并且用( λi , xi)作为目标特征对的近似，停止迭代。否则 只计算 Hm 的全部特征值 λi ，转第(4)步。
(4)选择 k 个 Ritz 值 λi 作为目标特征值的近似。每个处理 器进程用式式(1)式和式式(2)式计算 k/P 个精化向量及残量范
基于传统的重启动 Arnoldi 并行算法，本文提出用任务 图模型并行计算精化向量，避免了用数据并行模型使得通信 次数增加的问题。，从而使得本文提出的精化重启动 Arnoldi 并行算法具有较好的并行效率和可扩展性。
2 精化重启动 Arnoldi 算法
2.1 基本 Arnoldi 过程 给定 n 阶方阵 A 和单位初始向量 v1，基本 Arnoldi 过程
第 36 卷 第 6 期 Vol.36 No.6 ·博士论文·
计算机工程 Computer Engineering
文章编号：1000—3428(2010)06—0012—03
文献标识码：A
2010 年 3 月 March 2010
中图分类号：TP301.6
求解大规模矩阵特征问题的并行算法研究
赵 韬，迟学斌，陆忠华，赵永华
对于精化重启动 Arnoldi 算法，注意到每个处理器进程
本地都有矩阵 Hm 及其特征值。设需要求解 k 个特征对。如果 将计算一个精化向量作为一个任务，那么每个处理器进程都
可以根据式式(1)式和式式(2)式独立完成 k/P 个任务，而不需
要与其他处理器进程通信。如果( λi , ui)满足收敛条件，则停 止迭代。这时每个处理器进程本地拥有全部 Ritz 值和部分精
化向量。如果不收敛，那么使用精化向量和全局归约函数构
造 v1 进行重启动。精化重启动 Arnoldi 并行算法描述如下。 算法 2：
(1)设置目标特征对的个数为 k，基本 Arnoldi 过程的迭代
步数为 m，允许误差为 tol；给出第 1 次重启动的单位初始向
量 v1。 (2)用算法 1 并行计算 m 步基本 Arnoldi 过程。每个处理
目前，传统的重启动 Arnoldi 并行算法使用数据并行模 型[3-4]并行计算 Krylov 子空间正交基。这种并行算法是将矩 阵和向量按行分配到各个处理器进程上。每个处理器进程首 先只对本地矩阵和向量进行计算，如果需要其他处理器进程 上的数据，则再与其他处理器进程通信获取所需要的数据。 当所需计算的数据较少并且处理器进程数增多时，应用数据 并行模型的并行算法的通信次数会显著增加，降低了算法的 并行效率和可扩展性。
用矩阵形式可以表示为 AVm = Vm Hm + hm+1,mvm+1em*
或者
—12—
AVm = Vm+1H m
其 中 ， m<<n, Vm ∈ Cn×m 是 列 正 交 矩 阵 ； Hm ∈ Cm×m 是 上
Hessenberg 矩 阵 ； 且 VmHvm+1 = 0 , ， Vm+1=(Vm, vm+1) ； H m ∈ C(m+1)×m 的前 m 行是矩阵 Hm，其最后一行只有第 m 个元
向量是否收敛。如果 σ min 足够小，那么( λi , ui)即可作为目标
特征对的好的近似。，否则由精化向量构造 v1 进行重启动。
3 精化重启动 Arnoldi 并行算法
基本 Arnoldi 过程的主要计算部分是 Gram-Schmidt 正交 化过程[4-5]。为了降低全局通信次数和调用 2 级 BLAS 以有效
化算法构造 K(v1, A)的正交基。 算法 1：
for j = 1,…, m
1
w
r j+1
=
func
(A,v
r j
)
,
zr
=||w
|| r 2
j+1 2
,
r=1,2,…,
P,
reduce: z = z1+ z2+…+ zP，, broadcast: z;
2 iter =stract】Based on the parallel algorithm for restarted Arnoldi method by using data parallel model, a parallel algorithm for refined restarted Arnoldi method is presented. To impair the negative impact on parallel performance due to weak scalability, Tthe parallel algorithm presented implements parallel computing refined vectors by using task graph model and reduces the communication cost. As a result, the algorithm achieves parallel computing efficiently. Numerical experiments on KD-50-I show that the parallel algorithm presented performs good scalability and efficiency. 【Key words】matrix eigenvalue; Arnoldi algorithm; parallel computing; refined vector
||
(A
−
λi I )u
||2
其中，K(v1, A)表示由 v1, A v1,…, Am-1 v1 构成的子空间，通常
称 K(v1, A)为 Krylov 子空间。而 v1, v2,…, vm 是 K(v1, A)的一
组正交基[2,4]。 定 理 [2] 设 zi 是 矩 阵 H m − λi I 的 属 于 最 小 奇 异 值
(中国科学院计算机网络信息中心超级计算中心，北京 100190)
摘 要：基于数据并行的重启动 Arnoldi 并行算法，基于使用数据并行模型的重启动 Arnoldi 并行算法，提出一个精化重启动 Arnoldi 并行 算法。为了降低弱扩展性对并行性能的负面影响，该算法使用任务图模型并行计算精化向量，减少处理器进程之间的通信次数，有效地实 现并行计算。在 KD-50-I 万亿次机上的测试结果表明，该算法具有较好的可扩展性和并行 效率。 关键词：矩阵特征值；Arnoldi 算法；并行计算；精化向量
3
sr=
(Vjr
)
H
w
r j+1
,
r=1,2,…,
P;
4 reduce: H(1:j, j)=H(1:j, j)+ s1+ s2+…+sP,