CN113704691B - 一种申威众核处理器的小规模对称矩阵并行三对角化方法 - Google Patents

Abstract

本发明公开了一种申威众核处理器的小规模对称矩阵并行三对角化方法，其步骤包括：1)对待处理的对称矩阵进行二维划分，得到多个子块；每个从核负责处理一个对应子块；2)对于该从核组中的每一列从核，第i列从核对其LDM空间中存储的列数据进行householder变换，并计算参数tau及向量v；将向量v存储到第i列从核的每一从核LDM空间中，并调用swblas提供的列规约函数计算对应列向量的2范数；3)该从核组产生用于对该对称矩阵更新的向量w；4)该从核组根据向量w和向量v对该对称矩阵进行秩2更新；5)对该对称矩阵的每一列重复进行步骤2)到步骤4)，完成该对称矩阵三对角化。

Description

一种申威众核处理器的小规模对称矩阵并行三对角化方法

技术领域

本申请涉及申威众核处理器的稠密线性代数并行计算领域，具体涉及一种国产申威众核处理器的对称矩阵并行三对角化方法。

背景技术

对称矩阵的特征值、特征向量计算是稠密线性代数领域非常重要的核心计算之一，在第一性原理计算、量子模拟、人工智能等领域使用广泛。由于对称三对角矩阵具有很多非常好的性质，实际计算对称矩阵的特征值时往往先将对称矩阵化成对称三对角阵，然后使用二分法、QR迭代、分治法等算法求解特征值。对称矩阵三对角化在整个特征值计算中耗时较长，是并行和优化的重点。

对称矩阵三对角化sytd2是稠密线性代数领域著名的开源库LAPACK中一个子函数。通过householder变换将对称矩阵化成三对角矩阵。其通过调用LAPACK函数larfg计算householder变换的参数tau和向量v，通过调用BLAS库函数symv(对称矩阵向量乘)、syr2(对称矩阵秩2更新)、axpy(向量更新)、dot(向量点积)来进行矩阵更新。国产申威众核处理器上提供了并行版BLAS库swblas，其针对申威众核处理器的架构特点进行了并行和优化。sytd2可通过调用swblas来实现并行和加速。然而，该并行方式加速sytd2效果非常差。因为该方式需要重复加载数据到从核的局部存储器中，产生了大量的访存开销。提高sytd2的性能，可大幅减少对称矩阵求特征值的时间，对实际应用有重要的意义。

由于sytd2的重要性，国内外很多专家学者对其开展了深入的研究。有些工作通过对算法进行改进(将其约化改成两步，第一步先约化成带状对称矩阵，然后再约化成三对角矩阵)来改善其多核性能，有些研究面向多核平台、GPU加速卡的并行算法及优化方法，均取得了不错的效果。而在国产申威众核处理器上尚未出现该算法的并行和优化方法，本发明主要是填补此项空白。

国产申威众核处理器是由我国自主研制的一款高性能处理器，最新型号为39000，其采用主从异构架构，每个处理器节点由6个核组组成，每个核组包括1个主核和64个从核，从核以8*8阵列排布。每个从核可发起一个线程。最高可由64个线程同时进行计算。每个从核拥有256KB的私有LDM(Local Data Memory，局部数据存储器)空间，可通过DMA(直接存储器访问，Direct Memory Access)将数据从主存搬到LDM中。另外从核间可通过RMA(远端内存访问，Remote Memory Access)来进行数据传输。

发明内容

本发明的目的是提供一种国产申威众核处理器上小规模对称矩阵三对角化的并行和优化方法，以解决上述现有技术运行时间过长的问题，减少对称矩阵三对角化在国产申威众核处理器上的运行时间。在实际应用中，需要计算特征值的矩阵往往比较小，所以本发明主要面向小规模的矩阵来设计，具体规模大小以能全部装入从核LDM空间为限。

为实现上述目的，本发明提供了如下方案：

本发明提供一种国产申威众核处理器上小规模对称矩阵三对角化的并行和优化方法。为了方便描述，首先给出对称矩阵三对角化的原始算法(SYTD2)：

for J＝1to N-1do

对当前列产生householder变换因子tau和向量v(larfg)

产生用于矩阵更新的向量w(symv、dot、axpy)

矩阵更新(syr2)

end

本发明提出的新方法是面向国产申威众核处理器从核组的，新方法相比原始算法更加复杂，需要进行数据加载、划分等。本发明所提出的方法只对矩阵进行一次加载，后面操作全部通过从核RMA进行数据交换，以完成所有计算任务。对称矩阵即是本文要处理的数据，后续所有操作均在矩阵数据上完成，具体步骤与原始算法类似，会增加一个数据划分的步骤，该步由所有从核协同完成，除了该步外，其他步骤均需由若干从核协同完成，具体内容如下：

(1)、数据划分。对于二维对称矩阵而言，并行数据划分通常有两种方式：一维数据划分和二维数据划分。一维数据划分会导致有些操作(如向量2范数计算)串行执行，从而性能较差。本发明采用二维数据划分，每个从核负责一个子块，见图3。初始时，由所有从核通过一个共同的DMA操作将矩阵数据从主存拷贝到各自的LDM空间，每个从核拥有划分后的一个子块。由于该操作面向的是对称矩阵，原始输入时只存储上三角或者下三角部分。在初始从内存将数据传输到从核LDM时，从核采用整块传输的方式获取对应子块，每个从核均拿到相应存储空间的数据。这样部分从核拿到的是无效数据，实际并不使用。这种方式DMA时数据更加规整，传输带宽利用更高，所以性能更高。

(2)、对当前列进行householder变换，并计算参数tau及向量v。由于矩阵采用二维划分，而该步需要使用当前矩阵的一列数据(其分别存储在一列从核的LDM空间中，该列从核称为当前从核列)，所以由该从核列协同完成。这种并行方式可有效减少该操作的计算时间。参数tau由larfg函数计算得到，向量v由当前列从核分别更新其所存储的部分得到。计算得到的向量v并不单独存储，而是与矩阵的当前列共享同一块空间，分别存在当前从核列的LDM中，以备后续计算使用。该步计算时，还需要调用swblas提供的列规约函数完成，以计算当前列向量的2范数。

(3)、产生用于对该对称矩阵更新的向量w。该步操作需要调用symv、dot、axpy来完成。symv主要完成对称矩阵和向量的乘法，对称矩阵的数据即步骤(1)中传输到每个从核LDM空间的数据信息，向量是步骤(2)中计算得到的向量v。symv操作由拥有数据的各个从核协同完成，但是向量v只有当前列从核拥有数据，其他从核计算需要的向量子块需要通过RMA传输。实际计算时，矩阵的对角线部分所在从核只保存了矩阵的上/下三角的数据，需要补全后再进行乘法操作。而对于矩阵的其他部分，只有部分从核拥有有效数据。以矩阵为下三角为例，只有位于从核阵列左下的从核拥有有效数据，见图2，其余从核拥有的是无效数据，图2中以空白示意。对于位于左下的从核，在进行矩阵向量乘计算时不仅要计算当前子块的部分，还要计算其对称部分相应的矩阵向量乘结果，这样可以重用矩阵的数据，减少数据传输开销。为此本发明在位于左下的从核设置两块缓冲区，用于存储中间矩阵和向量相乘的计算结果。当其两个部分即当前块部分和其对称部分的数据均计算完成后，需要将对称部分的中间结果传到对应从核，即位于左下的从核需要将其计算的对称部分结果传给位于右上的从核，以方便进行symv结果的累加。累加工作通过调用swblas提供的行规约函数完成。

由于位于左下的从核需要计算两个部分，相应的向量v也需要两个部分，RMA时需要分别进行传输。另外位于对角线的从核只需要一部分向量v，单独进行传输，具体传输模式见图3。

symv计算完成后，再通过dot和axpy来更新向量w，dot和axpy均由一列从核完成。

(4)、对矩阵进行秩2更新，由syr2完成。更新时需要传输向量w和向量v，由对角线及其右上/左下的从核进行矩阵更新。位于对角线的从核只需要传输向量w和v的一个子块，位于左下/右上的从核需要更新当前拥有的子块，但是该操作需要向量v和w的两个部分，具体可有syr2算法推导而来。所以向量的传输方式与symv操作类似，见图3。

(5)、上述步骤中，所需使用的BLAS和LAPACK函数larfg、dot、axpy、symv、syr2，均不再直接调用库函数，而是实现了单从核版，单从核版即是串行实现，与原始算法的计算过程一致。

(6)、上述步骤中，矩阵只加载一次，其余计算中所需要的数据均通过RMA来完成。国产申威众核处理器提供了大量的接口来实现从核间数据交换，包括rma_get、rma_row_bcast等，上述步骤中的数据交换，均通过调用相应接口完成。

(7)、为了尽量减少sytd2函数运行时间，本发明还进行了性能调优，对步骤(3)中的symv、步骤(4)中的syr2进行了向量化。

本发明的技术方案为：

一种申威众核处理器的小规模对称矩阵并行三对角化方法，其步骤包括：

1)对待处理的对称矩阵进行二维划分，得到多个子块；从核组中的每个从核负责处理一个对应子块；该从核组中的每列从核的LDM空间中存储该对称矩阵中的一列数据；

2)对于该从核组中的每一列从核，该从核组中的第i列从核对其LDM空间中存储的列数据进行householder变换，并计算参数tau及向量v；其中参数tau由larfg函数计算得到，向量v由第i列从核分别更新其所存储的子块数据得到；将所得向量v存储到第i列从核的每一从核LDM空间中，以及调用swblas提供的列规约函数计算对应列向量的2范数；

3)该从核组产生用于对该对称矩阵更新的向量w；

4)该从核组根据向量w和向量v对该对称矩阵进行秩2更新；

5)对该对称矩阵的每一列重复进行步骤2)到步骤4)，完成对该对称矩阵三对角化。

进一步的，该从核组的所有从核通过一个共同的DMA操作将该对称矩阵的数据从主存拷贝到各自的LDM空间，每个从核拥有划分后的一个子块。

进一步的，从核采用整块传输的方式获取对应子块。

进一步的，步骤2)中，将参数tau通过RMA广播给该从核组的所有从核。

进一步的，产生所述向量w的方法为：首先通过RMA将向量v传输给每个参与计算的从核；然后所有参与计算的从核分别根据向量v进行symv计算，得到向量y；然后该从核组的每一列从核分别通过dot计算该列从核中各从核所得向量y的部分和并对其进行规约后调用axpy完成向量w的更新。

进一步的，对存储有有效数据的从核a设置两块缓冲区，分别用于存储从核a对其所存储子块与向量v的相乘结果a1、从核a对其所存储子块在该对称矩阵中的对称部分与向量v的相乘结果a2；然后将a1作为从核a对应的向量y，将a2传输给该对称部分对应的从核b作为从核b对应的向量y。

进一步的，对该对称矩阵进行秩2更新的方法为：首先通过RMA将向量w和向量v传输给每个参与进行秩2更新的从核，然后由对角线及存储有有效数据的从核对该对称矩阵进行秩2更新。

本发明与现有技术相比的有益效果：

设计了一套SYTD2的众核并行算法，适用于国产申威众核处理器。该算法大量减少了数据搬运时间，有效利用了国产申威众核处理器的计算和访存资源，大幅减少了该函数的实际运行时间，相比原始版本，平均加速比达到了5倍左右。

附图说明

图1为本发明方法整体流程示意图；

图2为矩阵任务划分示意图；

图3为symv操作前向量RMA传输示意图；其中(a)是对角子块所需数据传输示意图，(b)是本地子块计算所需数据传输示意图，(c)是对称子块计算所需数据传输示意图。

图4为本发明新方法的效果示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

步骤一、数据划分，见图2。将对称矩阵进行二维划分，每个从核负责其中一个子块，子块数据存放在从核LDM中；利用国产申威众核处理器提供的DMA操作将数据分别从主存加载到从核LDM空间中。

步骤二、对当前列进行householder变换，并计算参数tau及向量v。首先由当前从核列分别计算部分和，然后将结果规约到当前从核(当前处理列起始元素所在的从核)。由当前从核计算向量2范数，并调用larfg函数计算参数tau。然后利用该列从核分别更新向量v。向量v由对称矩阵的当前列更新而来，更新后仍然存储在对称矩阵的当前列，以便后续使用。故向量v的数据是存放在所有当前列从核的LDM空间中的。计算的参数tau还需要通过RMA广播给所有从核，以便后续使用。

步骤三、产生用于矩阵更新的向量w。首先通过RMA传输每个参与计算的从核需要的向量v，具体需要的子块由对称矩阵向量乘的算法决定。具体传输分为3类：位于对角线的从核为一类，见图3(a)，将上步中计算的向量v每个子块传给当前子块所在从核行对应的对角从核上；第二类为传输当前子块矩阵向量乘需要的向量v子块，见图3(b)，将向量v的每个子块传给列号与该子块所在从核行号相等的从核；第三类为传输对称子块矩阵向量乘需要的向量v子块，见图3(c)，将向量v的每个子块传给其所在从核同行的其他所有从核。然后由所有参与计算的从核分别计算当前子块矩阵和向量v子块相乘的结果作为中间结果，其中位于对角线的从核需要先补全数据再进行计算。以矩阵为下三角为例进行说明，对于位于左下的从核，在进行矩阵向量乘计算时不仅要计算当前子块的部分，还要计算其对称部分相应的矩阵向量乘结果，这样可以重用矩阵的数据，减少数据传输开销。为此本发明在位于左下的从核设置两块缓冲区，用于存储中间矩阵(即子块)和向量相乘的计算结果。当其两个部分即当前块部分和其对称部分的数据均计算完成后，需要将对称部分的中间结果传到对应从核，即位于左下的从核需要将其计算的对称部分结果传给位于右上的从核，以方便进行symv结果的累加。累加工作通过调用swblas提供的行规约函数完成。

symv计算完成后，产生结果向量y，再通过dot和axpy来更新向量w，dot和axpy均由当前从核列完成。首先由该列从核分别计算向量y的部分和，由当前从核负责规约，并将计算得到的向量点积结果传给当前从核列的所有从核，向量w通过向量y和向量v，以及计算得到的向量点积结果来得到，该操作通过调用axpy完成。最终向量w也存放在当前从核列，每个从核存一个子部分。

步骤四：对矩阵进行秩2更新。首先通过RMA传输每个从核需要的向量w和向量v，具体传输分为3类：位于对角线的从核为一类，见图3(a)；第二类为传输当前子块矩阵更新需要的向量v、w子块，见图3(b)；第三类为传输当前子块矩阵向量乘需要的转置部分向量v、w子块，见图3(c)。然后由对角线及其存储有效数据的从核进行矩阵更新。

步骤五：对矩阵的每一列重复进行步骤二到步骤四，即可完成对称矩阵三对角化，最终对称矩阵只有三条对角线。

采用以上介绍的技术，用国产申威39000众核处理器进行验证。由于LDM所限，只有1024规模及以下的矩阵可以使用本发明来进行计算，测试数据见图4，(以单精度为例，实际本发明适用于所有其他精度计算)，从图中可以看出，测试选取的4个规模计算时间相比原始版本都有了明显的降低，平均加速比达到5倍左右。

以上所述的实施例仅是对本发明的优选方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。

Claims (7)

1.一种申威众核处理器的小规模对称矩阵并行三对角化方法，其步骤包括：

1)对待处理的对称矩阵进行二维划分，得到多个子块；申威众核处理器的从核组中的每个从核负责处理一个对应子块；该从核组中的每列从核的LDM空间中存储该对称矩阵中的一列数据；

2)对于该从核组中的每一列从核，该从核组中的第i列从核对其LDM空间中存储的列数据进行householder变换，并计算参数tau及向量v；其中参数tau由larfg函数计算得到，向量v由第i列从核分别更新其所存储的子块数据得到；将所得向量v存储到第i列从核的每一从核LDM空间中，并调用swblas提供的列规约函数计算对应列向量的2范数；

3)该从核组产生用于对该对称矩阵更新的向量w；

4)该从核组根据向量w和向量v对该对称矩阵进行秩2更新；

5)对该对称矩阵的每一列重复进行步骤2)到步骤4)，完成对该对称矩阵三对角化。

2.如权利要求1所述的方法，其特征在于，该从核组的所有从核通过一个共同的DMA操作将该对称矩阵的数据从主存拷贝到各自的LDM空间，每个从核拥有划分后的一个子块。

3.如权利要求1所述的方法，其特征在于，从核采用整块传输的方式获取对应子块。

4.如权利要求1或2或3所述的方法，其特征在于，步骤2)中，将参数tau通过RMA广播给该从核组的所有从核。

5.如权利要求1所述的方法，其特征在于，产生所述向量w的方法为：首先通过RMA将向量v传输给每个参与计算的从核；然后所有参与计算的从核分别根据向量v进行symv计算，得到向量y；然后该从核组的每一列从核分别通过dot计算该列从核中各从核所得向量y的部分和并对其进行规约后调用axpy完成向量w的更新。

6.如权利要求5所述的方法，其特征在于，对存储有有效数据的从核a设置两块缓冲区，分别用于存储从核a对其所存储子块与向量v的相乘结果a1、从核a对其所存储子块在该对称矩阵中的对称部分与向量v的相乘结果a2；然后将a1作为从核a对应的向量y，将a2传输给该对称部分对应的从核b作为从核b对应的向量y。

7.如权利要求5所述的方法，其特征在于，对该对称矩阵进行秩2更新的方法为：首先通过RMA将向量w和向量v传输给每个参与进行秩2更新的从核，然后由对角线及存储有有效数据的从核对该对称矩阵进行秩2更新。

Images