稠密对称特征值问题存在于科学计算和工程模拟中的许多领域，例如计算材料和计算化学等科学计算应用[1-3].规约稠密对称矩阵为三对角形式是求解对称特征问题的关键一步，且占据了问题求解的大部分时间.因此对稠密对称矩阵三对角化算法进行加速计算是非常必要的.虽然在LAPACK(linear algebra package)[4]，ScaLAPACK(scalable linear algebra package)[5]等一些传统的标准库中都可以获得对称矩阵三对角化算法，但硬件发展促使我们重新设计算法以更有效地使用底层硬件.本文基于GPU和CPU的混合架构重新设计了对称矩阵三对角化的并行算法的实现方法，使其可以有效地利用GPU的并行计算能力，降低特征值求解的计算时间.

稠密对称矩阵三对角化一般分为2种方法：1阶段法和2阶段法.标准的1阶段法是直接将稠密对称矩阵通过双边的Householder正交变换转化为三对角矩阵[4].而2阶段法[6-7]包括两步,第1步通过一系列的Householder变换将矩阵规约成一个带状(band)矩阵，这一步运用了许多BLAS(basic linear algebra subroutine)-3(矩阵-矩阵)操作，因此取得了较好的效果[7-9]，第2步将带状矩阵转化为三对角矩阵[10].在2阶段法中，由于第1步使用了更多的BLAS-3操作，通常2阶段法要比1阶段法更快.然而，由于2阶段法中第2步的计算不能充分发挥GPU的性能，并且若要求解原始问题特征向量，2阶段法比1阶段法多需要1个额外的转化计算，因此在整体上所需的时间要比1阶段法更多[11].因此，在本文中，我们采用1阶段法.

在1阶段法的三对角化块算法中，主要分为2部分：条块内规约与条块外尾部矩阵的修正.主要的计算量主要集中在条块内规约时的对称矩阵向量积和条块外尾部矩阵的修正.已经有一些工作利用GPU对三对角化中这2部分加速计算[11-12]，文献[11-12]采用了仅对尾部对称矩阵的上或下三角部分进行修正的方法，该方法降低了修正过程中的计算量，但使得修正后的矩阵失去了完整性，导致了条块内规约过程中对称矩阵向量积计算变为2个不规则的三角矩阵向量积计算，这种不规则的计算不能充分利用GPU的并行计算能力.针对GPU集群的架构特点，本文给出的算法改进包括4个方面：

1) 给出的三对角化算法是对整个矩阵进行操作，这避免了条块内规约过程中不规则矩阵向量乘运算,但这也带来了一些冗余计算：在条块内修正为规约部分时，不仅要修正列条块内为规约部分，也要修正行条块内未规约部分.但从整体上减少了算法的用时，增加了性能.

2) 增加计算密集度.在条块内规约中将2个小的长条矩阵-向量积合并成1个更大的矩阵-向量积.并且，在条块外尾部矩阵的修正中，将调用2次矩阵-矩阵乘的秩-b修正合并成1个调用1次矩阵-矩阵乘秩-2b修正.加大计算密集度能更加充分发挥GPU的计算能力，增加GPU的利用率，从而提升算法的性能.

3) 利用多个CUDA流同时启动多个CUDA内核，使得在三对角化算法中独立的CUDA操作可以在不同的CUDA流中并发执行.另外，在基于MPI+CUDA的并行算法中，利用了CPU与GPU之间的异步数据传输操作，使得在不同CUDA流中的数据传输操作与GPU内的核函数同时执行，隐藏了数据传输的时间，使得并行算法的性能进一步提高.

4) 采用完全并行的点-点数据通信方式.在并行算法中，由在处理器列上的向量获得对应的在处理器行上的向量操作，若使用全局通信，则需要将属于在处理器列上的各个列处理的部分向量进行全局求和操作，这是一个全局的非并行的通信操作.本文给出的处理器行列之间的向量转置算法是一个完全并行的点-点的并行通信形式：处理器列上的各个处理器并行地将属于本处理器的部分向量发送到对应的处理器行上的处理器，处理器行上的各个处理器并行地接收发送来的数据.这有效地改善了三对角化算法实现的通信性能.

本文三对角化算法的实现是在我们开发的PSEPS(parallel symmetric eigenproblem solver)特征求解器软件包的基础上完成的.我们从多个方面优化了三对角化算法实现，使其主要计算部分并行地在多个GPU上进行加速计算.计算部分通过调用NVIDIA的cuBLAS(CUDA basic linear algebra subroutine)库完成，从而有效地保证了多GPU实现的性能.在中国科学院超级计算机系统“元”上，使用Nvidia Tesla K20 GPGPU卡(峰值双精度浮点性能1.17 TFLOPS)测试了算法的实现性能，取得了较好的加速效果与可扩展性.在对30 000×30 000阶矩阵利用GPU进行加速计算时，文献[8]中利用4个NDIVIA Tesla C2050 GPUs(峰值双精度浮点性能515.2GFLOPS)加速时，性能约为115GFLOPS，约达到峰值性能的22.32%.我们利用4个Nvidia Tesla K20 GPUs时性能约为376.60 GFLOPS，约达到峰值性能的32.19%.

1 三对角化串行块算法

定理1[13].给定向量a=(a1,a2,…,an)∈n，可以得到向量u∈n和常量β，使得:

其中，a1是a的第1个元素，e1是第1个单位向量.事实上和常量β=是满足定理1的数量对.我们称定理1中的(I-β uuT)为Householder变换；u称为对应于a的Householder向量.为了计算Householder变换，需求解(u,β,μ)，记为：(u,β,μ)=House(a).

对于给定的实对称矩阵A∈n×n，可以通过对其执行n-2次双边Householder正交变换{H1,H2,…,Hn-2}，使得:

T=Hn-2…H2H1AH1H2…Hn-2

为三对角对称矩阵.

将矩阵其中a11是常量)进行一次两端Householder正交变换:

(I-β uuT) A(I-β uuT)=A-uwT-wuT，

其中该变换是一个对A进行秩-2修正的算法，主要的计算是BLAS-1(向量-向量)操作，而GPU适合计算密集型任务，BLAS-1计算并不能充分利用GPU的计算性能.因此需重新设计算法使修正过程包含较多的BLAS-3操作.常用的方法是块算法(blocked algorithm)，积累多个BLAS-1操作为BLAS-3操作.本文假设块的大小为b，将矩阵A划分成一系列的宽度为b的条块，则三对角化的块算法可以分成如图1所示的2部分：图1(a)是条块内规约，该操作对条块内的每一列进行规约，产生的向量u和w分别保存在宽度为b的矩阵U和W中；图1(b)是尾部矩阵A22的修正，该操作在条块规约完成后进行，利用在条块规约中产生的矩阵U和W，对尾部矩阵A22进行修正，A22=A22-UWT-WUT.这将多个BLAS-1积累为BLAS-3操作.之后在对尾部矩阵A22进行相同的操作，直到不存在条块为止.

块算法将尾部矩阵的秩-2修正聚合为秩-b修正.A22的延迟修正使向量v的生成变得复杂，因为存储在矩阵A22中的元素在条块规约过程中没有及时更新.因此，当求解w时，v的计算需要替换为

v=β(A22-UWT-WUT) u，

其中,UWTu与WUTu对应着尾部矩阵的修正.进一步：

v=β A22u-β(UWT+WUT) u=

β A22u-β MNTu，

其中,M=(U W),N=(W U)∈n×2b.这将原本需要4次矩阵向量乘的运算合并，只需要2次即可完成.

同样，在尾部矩阵A22的修正中：

A22=A22-UWT-WUT=A22-(WUT+UWT)=

这将需要2次矩阵乘法的秩-b修正合并成了只需要1次矩阵乘法的秩-2b修正.在实现中，我们可以将向量w和u交叉存储在矩阵M和N中：

M=(u1,w1,u2,w2,…,ub,wb)，
N=(w1,u1,w2,u2,…,wb,ub).

将这些运算合并成一个更大的计算，这使得计算更加密集，更适合GPU的架构.密集计算可以更加充分地利用GPU的计算能力，从而提升了三对角化算法的性能.

算法1描述了利用1个GPU进行加速计算的对称矩阵三对角化块算法，其中Step1～Step7表示7个主要的部分，本文将对这7部分中的主要计算部分利用GPU进行加速.为了更好地发挥GPU的计算性能，我们的算法将对整个矩阵A进行操作，而不是仅对A的下三角或者上三角部分操作，这使得数据规整，本需要2次矩阵向量积的Step2计算现在仅需要1次矩阵向量积运算，但这同时也会增加Step6中更新条块内未规约部分的冗余计算：不仅要更新列条块中未规约部分另外也需要更新行条块内的未规约部分但相比于Step2而言，这部分冗余计算所用的时间较少.因此，从整体上来说可以提升三对角化算法的性能.算法1产生的Householder向量u存储在矩阵A的下三角部分，三对角矩阵的主对角元素存储在矩阵A的主对角位置，常量β存储在向量φ中，次对角元素μ存储在向量z中.算法1中，变量前边加“d”表示在GPU端对应的变量.

算法1. 利用GPU加速计算的三对角化串行块算法.

划分矩阵,其中dATL是0×0阶矩阵，置i=1;

do until dABR是0×0阶矩阵{

确定列条块的大小：bc=min{b,dABR的列宽}；

重新划分矩阵,其中dA01，dA11，dA21的列数为bc:

置(其中是0×0阶矩阵):

不是b×b阶do {

/*条块内规约*/

重新划分，其中是常量:

/*将u从GPU传到CPU*/

Synchronize(stream1);

Step1. 计算(u,β,μ)=House(u):

β→φi;μ→zi;i=i+1;

(stream1) u→du;

/*将向量u传到GPU*/

/*将du保存在dA中*/

Step2. 计算矩阵向量积：

Step3. 修正矩阵向量积:

Step3.1. (stream2) dvt=dNTdu;

Synchronize(stream1)；

Step3.2. (stream2) dv=dv-β dMdvt;

Step4. 计算修正向量:

Synchronize(stream2);

Step5. 将du和dw交叉保存到矩阵dM和dN中：

Step5.1. (stream1) dM=[dM|du|*];

Step5.2. (stream2) dM=[dM|*|dw];

Step5.3. (stream3) dN=[dN|dw|*];

Step5.4. (stream4) dN=[dN|*|du];

Step6. 更新条块内未规约部分：

Step6.1. 更新列条块：

Step6.2. 更新行条块：

Synchronize(stream1&stream2&stream3

& stream4&stream5&stream6);

对下列矩阵形式继续上述过程：

}end while

Step7. 对尾部矩阵dA22执行秩-2b修正：

(stream1) dA22=dA22-dNdMT;

对下列矩阵形式继续上述过程:

}end do

在算法1内的条块内规约过程中，可以发现Step2和Step3.1以及Step5.1～Step5.4，Step6.1～Step6.2这些操作互不依赖，每个操作可以看成是1个独立的任务，这些任务可以并发执行.在利用GPU进行加速计算的三对角化算法中，可以将每一个CUDA流视为GPU上的1个任务，这样通过启动多个CUDA流，使得这些独立的任务在不同的CUDA流中并发执行，使时间进一步压缩，从而提升三对角化算法的性能.注意在利用多个CUDA流并发启动的三对角化算法中，在有数据依赖的位置需要进行同步操作，以保证算法的正确性，如图2所示:

2 基于MPI+CUDA的三对角化并行块算法

2.1 数据在处理器上的分布方式

MPI并行算法的效率较大地依赖数据在处理器间的移动代价,这个代价又与数据在处理器上的初始映射(分布)和处理器间通信所使用的算法有关,并且数据在处理器的分布方式极大地影响着算法的负载平衡和通信特性.在稠密线性代数中，当前普遍使用的是2维块循环分布.2维块循环具有支持可拓展的高性能并行稠密线性代数编码的特性[14-16]，并被HPF(high performance fortran handbook)[17]，ScaLAPCK[5]和其他并行稠密线性代数库选为首要的数据分布.因此在基于MPI+CUDA的并行三对角化块算法中，我们采用了2维块循环的数据分布方法.将稠密对称矩阵A采用2维块循环方式分布到 p=pr×pc处理器网格中，存储块的大小为b.这里，假设矩阵A的元素A(1,1)映射到处理器p(0,0)中的本地矩阵元素为Aloc(1,1)，则全局元素A(i,j)映射到处理器

中的本地矩阵元素为Aloc(iloc,jloc)，其中

iloc=b+mod (i-1,b×pr)+1，
jloc=b+mod (j-1,b×pc)+1.

2.2 MPI通信方式

在2维通信域中，处理器间的通信通常发生在处理器网格中同行处理器之间或同列处理器之间.在行、列映射方式下矩阵的整行或整列被播送.基于MPI+CUDA并行三对角化算法中的通信主要操作包括行播送、列播送、行聚集求和、转置操作.

1) 行播送.2维通信域中处理器行内某处理器将数据广播到该处理器所在行的所有处理器，在算法2中此操作用Rbroadcast表示.

2) 列播送.2维通信域中处理器列内某处理器将数据广播到该处理器所在列的所有处理器，在算法2中用Cbroadcast表示.

3) 行聚集求和.在处理器行上对数据进行聚集求和，在算法2中用Rsum表示.

4) 转置操作.将处理器行/列上的向量转置到对应的处理器列/行上的向量.在算法2中用Transpose表示.

2.3 基于MPI+CUDA的并行块算法

基于串行的利用GPU加速计算的的三对角化块算法1，算法2给出了基于MPI+CUDA的并行三对角化块算法在GPU集群上的实现.在算法2中，(GC-数字)表示数据从GPU传到CPU，(CG-数字)表示数据从CPU传到GPU.(CC-数字)表示CPU之间的通信操作.(G-数字)表示在GPU端的操作.算法2中的变量前边加“d”表示在GPU端对应的变量.

在基于MPI+CUDA的并行算法2中，除了将独立的CUDA操作在不同的CUDA流中进行运算，还利用了CPU与GPU之间的异步数据传输操作，使得在不同CUDA流中的数据传输操和GPU内的核函数同时执行： (G-3)GPU内的拷贝操作与(GC-3)数据传输操作，(CG-4)数据传输操作与(G-4) GPU内的计算操作.这可以隐藏CPU与GPU之间数据传输的时间，使得基于MPI+CUDA的并行算法的性能进一步提高.

这里假设每个处理器拥有矩阵A的行号集合为RS，列号集合为CS.并且在三对角化之前已经将部分矩阵A(RS,CS)传到了对应的GPU中：A(RS,CS)→dA(RS,CS).为了使得多个GPU并行加速计算，这里采取的策略是CPU与GPU一一对应，因此GPU形成了一个与CPU一样处理器网格.这里采取条块的大小与存储块的大小相等，共要进行bn=n/b次条块规约.

算法2. 基于MPI+CUDA的并行三对角化块算法.

/*算法中的dot表示向量的点乘运算，scal表示向量的数乘运算，copy表示向量的拷贝操作，gemv表示矩阵-向量乘运算，axpy表示向量的加法运算，ger表示矩阵的一阶修正运算，gemm表示矩阵-矩阵乘运算*/

for each bt∈[1,bn]

as=bt×b;/*条块内第一列*/

ae=min(as+b-1,n-2);/*最后一列*/

for each k∈[as,ae] /*条块内规约*/

Step1. 计算Householder向量及其范数:

if k∈CS then CS=CS\(k) endif

if k∈RS then

Step1.1. (GC-1) (stream1) dA (k,CS)→uc;

Synchronize(stream1);

Step1.2. (CC-1) Cbroadcast(uc);

Step1.3. (CG-2) (stream1) uc→duc;

endif

/*dot*/

Synchronize(stream1);

Step1.5. (CC-2) t=Rsum(g);

if j+1∈CS then /*行播送*/

Step1.6. (CC-3) Rbroadcast(u(k+1));

endif

Step1.7. 计算次对角元素μ和范数β

Step1.8. (G-2) (stream1) duc=β duc;

/*scal*/

if k∈RS then

Synchronize(stream1);

Step1.9. (G-3) (stream2) duc→dA(k,CS);

/*copy*/

endif

if k∈RS then RS=RS\(k) endif

Step1.10. (GC-3) (stream1) duc→uc;

Synchronize(stream1);

Step1.11. (CC-4) ur=Transpose(uc);

if k+1∈CS then

Step1.12. (CC-5) Rbroadcast (ur);

Step1.13. (CG-4) (stream2) ur→dur;

endif

Step2.计算矩阵向量积(gemv):

(G-4) (stream1) dxr=β dA(RS,CS) duc;

Step3. 修正矩阵向量积:

Step3.1. (G-5) 计算部分向量dxr(gemv×2)

Synchronize(stream1&stream2);

(G-5.2) (stream0) dxr-=dMrdxt;

Step3.2. (GC-5) (stream0) dxr→xr;

Synchronize(stream0);

Step3.3. (CC-6) wr=Rsum(xr);

Step3.4. (CG-6) (stream2) wr→dwr;

Step3.5. (CC-7) wc=Transpose(wr);

if k∈RS then

Step3.6. (CC-8) Cbroadcast(wc);

Step3.7. (CG-7) (stream1) wc→dwc;

endif

Step4. 修正矩阵向量积:

/*dot*/

Synchronize(stream1);

Step4.2. (CC-9) σ=Rsum(s);

Step4.3. (G-7) 计算dwc, dwr(axpy×2);

Synchronize (stream1&stream2);

Step5.(G-8) 将向量dwr，dwc，dur，duc交叉保存到矩阵dMr,dNc的相应位置(copy×4):

(G-8.1) (stream1) dMr=[dMr|dur|*];

(G-8.2) (stream2) dMr=[dMr|*|dwr];

(G-8.3) (stream3) dNc=[dNc|dwc|*];

(G-8.4) (stream4) dNc=[dNc|*|duc];

Step6. (G-9)更新条块内未规约部分:

(G-9.1) (stream5) /*ger×2*/

(G-9.2) (stream6) /*ger×2*/

Synchronize(stream1&stream2&stream3&

stream4&stream5&stream6);

endfor

Step7. 对尾部矩阵执行秩-2b修正(gemm):

Synchronize(stream0);

endfor

算法2中CPU之间的通信包括如图3所示9次操作：1)2次Cbroadcast.(CC-1)和(CC-8)，分别需要lb pr次消息启动，传递的数据长度为(n-k)/pc;2)2次Rbroadcast.(CC-3)和(CC-5)，分别需

要lb pc次消息启动，(CC-3)传递的消息长度为1，(CC-5)传递的消息长度为(n-k)/pr;3)3次Rsum.(CC-2)，(CC-6)，(CC-9)，分别需要lb pc次消息启动，(CC-2)和(CC-9)传递的消息长度为1，(CC-6)传递的消息长度为(n-k)/pr;4)2次Transpose转置操作(CC-4)和(CC-7).该转置操作是一个点对点的并行通信方式，以(CC-4)为例进行分析，因为通过(CC-1)的列播送，则uc所在的处理列的每个列处理器已经包含了其完整拷贝，这样，处理器列中的各个列处理器可以并行地操作.因为uc采用的也是2维块循环分布，所以首先定位到属于这个列处理器的部分向量，然后将其发送到对应的处理器行上的行处理器，所有的列处理器并行的发送，行处理器并行地接收.

图4描述了将处理器列上的向量vc转置到对应的处理器行上的向量vr的过程，每一个列处理器并行的将部分向量发送到对应的行处理器上，行处理器并行地接收部分向量：pse,0→pre,0，pse,1→pre,1，pse,2→pre,0，pse,3→pre,3.将处理器行上的向量转置到对应的处理器列上的操作类似，这里不再讨论.

CPU+GPU的异构编程需要CPU与GPU之间的数据传输，要想利用GPU进行加速计算，需要将计算数据上传到GPU端.因为使用的多个GPU分布在多个节点上，则GPU间的通信需要在CPU中通过MPI进行，首先将要通信的数据从GPU端传回到CPU端，之后通过MPI通信，通信完成后，再将数据传回到相应GPU中.算法2中CPU与GPU的数据传输操作如图2所示，共有7次.

3 实验结果与分析

因为Frank矩阵

A=(aij),aij=n-max(i,j)-1,
i,j=1,2,…,n

具有解析特征值：

所以为了验证算法实现的准确性和测试算法的性能，我们使用了Frank矩阵作为测试数据.

我们在中国科学院计算机网络中心超级计算机系统“元”上实现了本文描述的算法，并在GPU队列上进行测试.每台GPGPU计算节点配置2块Nvidia Tesla K20 GPGPU卡(峰值双精度浮点性能(板卡)1.17 TFLOPS)，2颗Intel E5-2680 V2(Ivy Bridge|10C|2.8 GHz)处理器.使用的cuBLAS版本是6.0.37，使用的MPI是openmpi-1.8.4.在测试时，为了使得多个GPU并行加速计算，每个节点申请2个进程，1个进程对应1个GPU.测试的数据类型是双精度浮点数.

3.1 基于MPI+CUDA的三对角化并行块算法性能分析

表1展示了20 000×20 000阶矩阵三对角化块算法在CPU调用上MKL(math kernel library) 2013的BLAS库和通过多个CUDA流并发调用GPU的cuBLAS库使用不同进程数(1，2，4，8，16，32)时的时间对比.从表1中可以看出，首先利用GPU进行加速计算的三对角化算法取得了明显的加速效果.其次，随着使用GPU数目的增多，加速倍数在不断减小，这是因为随着随着使用GPU数目的增多，在单个GPU上负载减小，并且随着进程数目的增多，通信时间、CPU与GPU之间的数据时间等开销占比增大，因此计算加速效果不如负载大的好.总的来说，利用GPU对三对角化算法进行加速取得了较好的加速效果(这里b=64).

图5展示了当b=64时不同规模的矩阵使用不同GPU个数(1，2，4，8，16，32)时达到的性能.从图5中可以发现，在使用GPU的数目一样时，矩阵规模越大达到的性能越高，这是因为问题规模越小，GPU的负载越小，而GPU适合计算密集型任务，这样加速效果就不如规模大的好.并且规模越大，计算时间所占的比例越大，加速效果越明显，则达到的性能越高.从图5可以看出，基于MPI+CUDA的三对角化并行块算法取得了较好的性能.

图6展示了不同矩阵规模(10 000×10 000，15 000×10 000，20 000×20 000，25 000×25 000，30 000×30 000)使用不同个数的GPU进行加速计算时相对于使用1个进程/GPU时加速比曲线.

从图6中可以看出当矩阵规模一定时，随着使用GPU数目增加，加速比增加的幅度在不断降低，这是由于随着使用GPU个数增加，GPU与CPU之间的数据传输、通信等开销所占的比例增加导致的.当使用的GPU的数目一样时，较大规模矩阵的加速比大于较小规模矩阵的加速比，加速比随着矩阵规模的增大而增大，这是因为大规模的矩阵在GPU上的负载较大，比小规模矩阵更能够充分地利用GPU的计算能力.并且，随着使用GPU数目的增多，较大规模的加速比与较小规模的加速比相差越来越大，这是因为在单个GPU上的负载越大，取得的性能越好.总的来说，基于MPI+CUDA的三对角化并行块算法具有良好的可扩展性.

3.2 加大计算密集度的算法性能分析

不同的条块(矩阵存储块)的大小会影响三对角化算法的性能，当三对角化算法只启动1个CUDA流时，分别测试了如表2所示的当b分别取值32，64，128时的不同规模的矩阵(10 000×10 000，20 000×20 000，30 000×30 000)利用不同GPU个数(2，4，8，16)加速时的用时.可以看出，当b=64时用时最短，因此在只启动1个CUDA流的三对角化算法的分析中，取b=64.

表3展示了在三对角化算法中启动1个CUDA流时20 000×20 000阶矩阵使用1个GPU加速时各个计算部分在GPU上的用时.其中，算法2中(G-4)gemv计算条块内规约的矩阵向量积v=Au的计算量约为4n3/3，占据了在GPU上计算总时间的约86.62%，因此，采用完整的矩阵进行运算提升这部分的性能，将会提升整个程序的性能.

表4展示20 000×20 000阶矩阵的使用不同GPU个数(4，8，16，32)时分别采用对上三角和下三角直接调用2次cublasDgemv与使用完整矩阵只调用1次cublasDgemv的时间对比，可以发现采用完整矩阵的方式这部分时间减少了约一半，又因为这部分占了三对角计算时间的大部分，这部分性能的提升将会大幅度提升整个算法的性能.但采用完整矩阵的方式将增加(G-9)ger的计算量，表5展示了增加前后的时间对比，虽然这部分的用时增加了，但是相对于(G-4)gemv时间的减少这部分的增加的时间非常少，因此整个程序的性能相对提升了.

对条块内第(G-5)gemv步进行优化，将原来需要4次矩阵向量积的运算合并成只需要2次矩阵向量积的运算.表6展示了20 000×20 000阶矩阵的使用不同GPU个数(1，2，4，8，16，32)时合并前后的时间对比，从表6中可以发现，将矩阵-向量运算合并后，这部分的时间减少了大约42%.这是因为将小规模的矩阵向量乘运算合并成较大规模的矩阵向量乘运算后，更加充分地利用了GPU的资源，因此程序的性能得到了进一步地提升.

(G-10)gemm条块外尾部矩阵的秩-2b修正计算量约为4n3/3，虽然这部分的计算量与(G-4)gemv的计算量一样，但占据的时间仅为7.27%(20 000×20 000阶矩阵使用1个GPU加速时)，这2部分时间占比相差非常大是因为BLAS-2与BLAS-3对GPU的利用率不同导致的.在(G-10)gemm中我们将需要2次gemm的秩-b修正合并成调用1次gemm的秩-2b修正，表7展示了20 000×20 000阶矩阵使用不同GPU个数(1，2，4，8，16，32)时优化前后的时间对比，从表7可以发现，将秩-b修正合并为秩-2b修正后，此部分的时间约减少了23%.合并计算后使得计算更加密集，更加充分地利用了GPU的计算能力，提高了GPU的利用率，因此获得了更高的性能.

图7展示了在三对角化算法中启动1个CUDA流时20 000×20 000阶矩阵使用不同GPU数目(2，4，8，16，32)加速时三对角化算法各个部分所用占用的时间(compu_gpu表示在GPU上的计算时间；memcpy表示CPU与GPU之间的数据传输时间；ctime表示CPU之间的通信时间；others表示其他时间).从图7中可以看出,在整个三对角化的过程中，计算时间占据了大部分，因此利用GPU进行加速计算能够明显提升三对角化算法的性能.并且，随着GPU的数目增多，计算时间compu_gpu逐渐递减.同时，通信时间ctime随着处理器数目(从2个开始)的增加也在逐渐递减.总的来说，整个基于MPI+CUDA的三对角化算法的用时随着使用GPU数目的增多逐渐递减，整个算法表现出了良好的强扩展性.

3.2 利用多个CUDA流的算法性能分析

在通过多个CUDA流并发调用多个GPU内核的三对角化算法中，我们同样测试了条块(矩阵存储块)b取不同值(32，64，128)时对不同规模的矩阵(10 000×10 000，20 000×20 000，30 000×30 000)利用不同GPU个数(2，4，8，16)加速时的影响，如表8所示，当b=64时用时最短，因此在启动多个CUDA流的三对角化算法的分析中，取b=64.

表9展示了在基于MPI+CUDA的并行三对角化块算法中启动1个CUDA流与启动多个CUDA流时，20 000×20 000阶矩阵使用不同进程数(1，2，4，8，16，32)时的时间对比，可以发现启动多个CUDA流相比于只启动1个CUDA流用时更短，取得了较好的加速效果.另外随着进程数目的增加，加速效果越来越好，这是因为随着进程数目的增加，在单个GPU上的负载不断减小，GPU中可用的SM(streaming multiprocessor)等资源越多，多个流并发执行的可能性越大，因此规模越小，加速效果越好.

Table 9 Time Comparison Between Starting a CUDA Stream

and Multiple Streams of 20 000×20 000 Order Matrix

表9 20 000×20 000阶矩阵启动1个CUDA流与多个CUDA流的时间比较

表10展示了不同规模的矩阵(10 000×10 000，20 000×20 000，30 000×30 000，40 000×40 000，50 000×50 000)利用32个GPU进行加速计算时启动1个CUDA流与启动多个CUDA流时的时间比较，可以发现随着矩阵规模的增大，加速效果不断减小.这是因为，虽然从软件的角度来看，CUDA操作在不同的流中并发运行；而从硬件上来看，不一定总是如此.根据每个SM资源的可用性，完成不同的CUDA流可能仍然需要互相等待.矩阵规模的越大，GPU上的负载越大，而GPU中SM资源是有限的，规模越大，需要的资源越多，因此多个内核同时并发执行的可能性就越小.并且，多个CUDA流并发执行，所用的总时间主要去取决于用时最多的那个任务，利用多个CUDA流并发调用多个内核，并没有减少用时最长的任务所用的时间，在三对角化算法中，主要是条块内规约中的矩阵向量积(G-4)gemv操作，矩阵规模越大，这部分的计算时间占的比例越多，因此对于整个三对角化算法，利用多个CUDA流并发调用多个内核，问题规模越大加速效果反而不如规模小的明显.总的来说，在单个GPU上的负载越小，利用多个CUDA流进行并发操作提升的性能越明显.

4 结束语

本文描述了基于MPI+CUDA的稠密对称矩阵三对角化并行块算法在GPU集群上实现，主要的计算部分由GPU完成，而通信部分在CPU上通过MPI来完成，从而充分发挥CPU+GPU的异构平台的优势.为了适应GPU架构，我们重构了算法，在条块内规约过程中的对称矩阵向量相乘时使用了完整的矩阵，虽然带来了额外的计算，但比直接分别对上三角和下三角矩阵进行2次矩阵-向量相乘时时间明显减少.同时，在条块内规约中合并了一些矩阵-向量积运算，这部分的时间减少了约42%.而且，在条块外尾部矩阵的修正过程中，将矩阵的秩-b修正合并为秩-2b修正，此部分的时间约减少了23%.本文提出的三对角化算法通过加大计算密集度来发挥GPU的计算能力，增加GPU的利用率，从而提升了算法的性能.另外，在三对角化算法中，对在GPU中互不依赖的操作，利用了多个CUDA流同时启动多个内核，使这些独立的CUDA操作可以在不同的流中并发执行，使得时间进一步压缩.并且，在基于MPI+CUDA的并行三对角化算法中，利用了CPU与GPU之间的异步数据传输，使得在不同流中的数据传输与核函数同时执行，隐藏了数据传输的时间，进一步提升了并行算法的性能.在三对角化并行算法中，块的大小将会影响通信次数、负载平衡以及利用GPU计算的性能，在使用时需要权衡.总的来说，本文中的异构三对角化算法取得了较好的性能与加速效果，并且具有良好的可扩展性.

参考文献

[1]Kent P R C. Computational challenges of large-scale, long-time, first-principles molecular dynamics[J]. Journal of Physics: Conf Series, 2008, 125(1): No.012058

[2]Auckenthaler T, Blum V, Bungartz H J, et al. Parallel solution of partial symmetric eigenvalue problems from electronic structure calculations[J]. Parallel Computing, 2011, 37(12): 783-794

[3]Singh D J, Nordstrom L. Planewaves, Pseudopotentials, And The LAPW Method[M]. Berlin: Springer Science & Business Media, 2006

[4]Anderson E, Bai Z, Bischof C, et al. LAPACK Users’ Guide[M]. Philadelphia: SIAM, 1999

[5]Blackford L S, Choi J, Cleary A, et al. ScaLAPACK Users’ Guide[M]. Philadelphia: SIAM, 1997

[6]Bischof C, Van Loan C. The WY representation for products of Householder matrices[J]. SIAM Journal on Scientific and Statistical Computing, 1987, 8(1): s2-s13

[7]Bischof C H, Lang B, Xiaobai S. The SBR Toolbox—Software for successive band reduction[J]. ACM Transactions on Mathematical Software, 2000, 6(4): 602-616

[8]Gansterer W N, Kvasnicka D F, Ueberhuber C W, et al. Multi-sweep algorithms for the symmetric eigenproblem[C] //Proc of Int Conf on Vector and Parallel Processing. Berlin: Springer, 1998: 20-28

[9]Dongarra J J, Sorensen D C, Hammarling S J, et al. Block reduction of matrices to condensed forms for eigenvalue computations[J]. Journal of Computational and Applied Mathematics, 1989, 27(1/2): 215-227

[10]Haidar A, Ltaief H, Dongarra J, et al. Parallel reduction to condensed forms for symmetric eigenvalue problems using aggregated fine-grained and memory-aware kernels[C] //Proc of 2011 Int Conf for High Performance Computing, Networking, Storage and Analysis. New York: ACM, 2011: 1-11

[11]Yamazaki I, Dong T, Tomov S, et al. Tridiagonalization of a dense symmetric matrix on multiple GPUs and its application to symmetric eigenvalue problems[J]. Concurrency and Computation: Practice and Experience, 2014, 26(16): 2652-2666

[12]Yamazaki I, Dong T, Tomov S, et al. Tridiagonalization of a symmetric dense matrix on a GPU cluster[C] //Proc of 2013 IEEE Int Symp on Parallel & Distributed Processing, Workshops and Phd Forum. Piscataway, NJ: IEEE, 2013: 1070-1079

[13]Demmel J W. AppliedNumerical Linear Algebra[M]. Philadelphia: SIAM, 1997

[14]Lichtenstein W, Johnsson S L. Block-cyclic dense linear algebra[J]. SIAM Journal on Scientific Computing, 1993, 14(6): 1259-1288

[15]Bolen J, Davis A, Dazey B, et al. Massively parallel distributed computing: World’s first 281 gigaflop supercomputer[C/OL] //Proc of the Intel Supercomputer Users Group. 1995 [2019-06-01]. http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=254434BE9D0FC733C1928A4E92BCF8C1?doi=10.1.1.26.3686&rep=rep1&type=pdf]

[16]Bisseling R H, van de Vorst J G G. Parallel LU decomposition on a transputer network[G] //LNCS 384: Parallel Computing 1988. Berlin: Springer, 1988: 61-77

Efficient Implementation of Parallel Symmetric Matrix Tridiagonalization Algorithm on GPU Cluster

Liu Shifang1,2, Zhao Yonghua1, Yu Tianyu1,2, and Huang Rongfeng1,2

1(Computer Network Information Center, Chinese Academy of Sciences, Beijing 100190)

2(University of Chinese Academy of Sciences, Beijing 100190)

Abstract The symmetric matrix tridiagonalization is the key computational process for solving dense eigenproblems. This paper presents the implementation of the dense symmetric matrix tridiagonal hybrid parallel blocked algorithm based on MPI(message passing interface)+CUDA(compute unified device architecture) for GPU cluster. The parallel algorithm design uses a two-level parallel method of MPI cluster level and GPU level. In the MPI-level parallelism, the communication performance of the tridiagonal MPI parallel algorithm is improved by designing the global data communication between the row-column communication domains in the two-dimensional communication domain as a completely parallel point-to-point data communication method. Moreover, by improving the original matrix tridiagonalized MPI parallel algorithm, the use of irregular matrix-vector operation in GPU-level parallelism is avoided, and the parallel performance of this part is improved by about 1 time. What’s more, the small-grained computing existing in the GPU parallel is merged into a larger granularity calculation, which fully utilize the computing power of the GPU by increasing the computational intensity, thereby increasing the utilization of the GPU and improving the performance of the algorithm. In addition, multiple CUDA streams can be used to enable independent CUDA operations in the algorithm to be concurrently executed in different streams. Furthermore, in the parallel algorithm, the asynchronous data transmission between the CPU and the GPU is utilized, so that the data transmission and the kernel function in different streams are simultaneously executed, which hides the time of data transmission and improves the performance of the algorithm. On the supercomputer system Era of the Computer Network Center of the Chinese Academy of Sciences, each compute node is configured with 2 Nvidia Tesla K20 GPGPU cards and 2 Intel E5-2680 V2 processors, we tested the performance of the implementation of the tridiagonalization blocked algorithm with GPGPU cards. It has achieved better acceleration, performance and scalability.

Key words symmetric matrix tridiagonalization; MPI+CUDA; point-to-point data communication; computational intensity; CUDA streams; scalability

(liushifang@cnic.cn)

中图法分类号 TP301

DOI:10.7544/issn1000-1239.2020.20190731

收稿日期：2019-10-11；修回日期：2020-02-21

基金项目：国家重点研发计划项目(2017YFB0202202)；中国科学院战略性先导科技专项(C类)(XDC01040000)

This work was supported by National Key Research and Development Program of China (2017YFB0202202) and the Strategic Priority Research Program of Chinese Academy of Sciences (C) (XDC01040000).

通信作者：赵永华(yhzhao@sccas.cn)

[17]Koelbel C H, Loveman D B, Schreiber R S, et al. The High Performance Fortran Handbook[M]. Cambridge, MA: MIT Press, 1994Liu Shifang, born in 1994. PhD. Her main research interests include numerical parallel algorithm library implementation and optimization of parallel algorithm library under heterogeneous parallel machine.

Zhao Yonghua, born in 1966. PhD, professor, PhD supervisor. His main research interests include parallel algorithm and software, parallel programming model, spectral clustering data analysis and so on.

Yu Tianyu, born in 1985. PhD. His main research interests include numerical parallel algorithm and solver.

Huang Rongfeng, born in 1990. PhD. His main research interest is the efficient implementation of the basic math library in heterogeneous platforms.