CN113641956B - 面向SW26010-Pro处理器的1、2级BLAS函数库的高性能实现方法 - Google Patents

Abstract

本发明公开一种面向SW26010‑Pro处理器的1、2级BLAS函数库的高性能实现方法，包括：对问题进行任务划分，产生若干子问题，其中所述问题的结构包括向量、普通矩阵、对称矩阵或三角矩阵；若为向量、普通矩阵或对称矩阵时，将各子问题的运算分配给相应线程；若为三角矩阵时，将子问题对角部分的运算分给0号线程，非对角部分的运算分配给其他相应线程；拼接各线程的运算结果，获取所述问题的解。本发明实现了BLAS 1、2级函数的并行化，解决了线程间的数据依赖问题，并通过自适应调优机制，进一步提升了函数的性能。

Description

面向SW26010-Pro处理器的1、2级BLAS函数库的高性能实现 方法

技术领域

本发明涉及基础线性代数库BLAS(Basic Linear Algebra Subprograms)实现的领域，尤其涉及一种面向SW26010-Pro处理器的1、2级BLAS函数库的高性能实现方法。

背景技术

BLAS是一种基础线性代数子程序库，主要包含向量和矩阵的基本操作，是最基础和最重要的数学库之一，广泛应用于科学计算、气象预报、天体物理等领域。BLAS库是很多专业软件的核心，其中，BLAS 1、2级函数会被几乎所有有关矩阵运算的应用以及稠密线性代数算法软件包(如LAPACK，ScaLAPACK)反复调用多次。在数值矩阵分析、深度学习等方面的实践表明，BLAS 1、2级函数对提高应用的运算速度和充分发挥高性能计算机的性能具有至关重要的意义。

BLAS 1、2级函数实现向量-向量、矩阵-向量操作，共包含30多个函数,且包含单精度、双精度、复数单精度和复数双精度四种类型。BLAS 1、2级函数具有访存密集的特征，其性能受限于系统访存带宽，且函数数量较多，函数涉及的矩阵在内存中有多种数据排布方式。如何对数据进行合理的划分，充分利用高效的访存模式，提高数据重用率，是BLAS 1、2级函数库高性能实现的极大挑战。

在BLAS 1、2级函数的高性能实现方面，国内外已经有相当多的研究成果。李毅等面向多核龙芯3A实现了二级BLAS函数库(李毅,何颂颂,李恺.多核龙芯3A上二级BLAS库的优化[J].计算机系统应用,2011,20(1):163-167.)。伴随着GPU加速器的快速发展，BLAS 1、2级函数在GPU上的优化工作近年来也成为研究热点，Jian Yin等在Nvidia GPU上利用寄存器分块方法实现了并行GEMV(Jian Y,Hui Y,Xu W,et al.Highly parallel GEMV withregister blocking method on GPU architecture[J].Journal of VisualCommunication&Image Representation,2014,25(7):1566-1573.)，Weizhi Xu等在NvidiaGPU上实现了一个针对GEMV的性能调优框架，针对GEMV的输入规模选择最优的算法(W.Xuet al.,"Auto-Tuning GEMV on Many-Core GPU,"2012IEEE 18th InternationalConference on Parallel and Distributed Systems,2012,pp.30-36,doi:10.1109/ICPADS.2012.15.)。

SW26010-Pro是具有异构架构的众核处理器。在基于SW26010-Pro众核处理器架构的神威新一代超级计算机上，目前还没有部署定制化的高性能BLAS1、2级函数库，且现有的开源数学库在该平台上性能较低，无法对应用提供有效的性能支撑。因此急需设计并实现一种面向该众核平台的高性能BLAS 1、2级函数库，以充分利用申威众核处理器的访存带宽，且满足上层应用对申威众核平台上高性能BLAS 1、2级函数的迫切需求。

发明内容

本发明提供一种面向SW26010-Pro处理器的1、2级BLAS函数库的高性能实现方法，以满足在SW26010-Pro众核处理器上对BLAS 1、2级函数的需求，解决现有的开源数学库性能较低的问题。

一种面向SW26010-Pro处理器的1、2级BLAS函数库的高性能实现方法，其步骤包括：

1)对问题进行任务划分，产生若干子问题，其中所述问题的结构包括向量、普通矩阵、对称矩阵或三角矩阵；

2)若为向量、普通矩阵或对称矩阵时，将各子问题的运算分配给相应线程；若为三角矩阵时，将子问题对角部分的运算分给0号线程，非对角部分的运算分配给其他相应线程；

3)拼接各线程的运算结果，获取所述问题的解。

进一步地，通过以下策略产生子问题：

1)对于向量，每一向量段视作一子问题x_i′，其中i′为向量段编号，0≤i′≤k-1，k为子问题的数量；

2)对于普通矩阵，每个行块视作一子问题A_i，其中i+1是矩阵的行编号，0≤i≤k-1；

3)对于对称矩阵，每一列块视作一子问题A_j，其中j+1是矩阵的列编号，0≤j≤k-1；

4)对于三角矩阵，每一行块视作一子问题A_i。

进一步地，所述问题的结构为向量时，通过以下步骤获取所述问题的解：

1)将子问题x_i′分配给相应线程T_i；

2)线程T₀计算子问题x₀的解y₀；

3)利用公式y_i←α×x_i′+y_i，各线程T_i计算得到解y_i，其中α为第一权重值；

4)拼接解y_i，获取所述问题的解y。

进一步地，所述问题的结构为普通矩阵时，通过以下步骤获取所述问题的解：

1)将子问题A_i分配给线程T_i，其中0≤i≤k-1，k为子问题的数量；

2)基于向量x′与子问题A₀，线程T₀计算得到解y₀；

3)利用公式y_i←α×A_i×x′+β×y_i，各线程T_i计算得到解y_i，其中α为第一权重值，β为第二权重值；

4)拼接解y_i，获取所述问题的解y。

进一步地，所述问题的结构为对称矩阵时，迪过以下步骤获取所述问题的解：

1)对每一子问题A_j进行划分，获取对角子矩阵D_j及下三角子矩阵L_ij，并将划分后的子问题A_j分配给线程T_j；

2)对向量x′进行划分，获取若干子向量x′_j；

3)将对角子矩阵D_j的上三角部分用对应的下三角子矩阵L_ij补齐；

4)各线程T_j基于对角子矩阵D₀与子向量x′₀、或上三角部分中相应的下三角子矩阵L_i0与子向量x′_j，计算得到解y₀或解L_0j；各线程T_j基于下三角子矩阵L_(j+1)j与子向量x′_j，计算得到相应的解y_(j+1)j；

5)对于对角子矩阵、下三角子矩阵及下三角子矩阵的对称部分，各线程T_j分别利用公式y_j←D_j×x′_j+y_j、y_i←L_ij×x′_j+y_i及y_j←L_ij×x′_i+y_j进行迭代求解，并拼接相应的子解，得到所述问题的解y。

进一步地，所述问题的结构为三角矩阵时，，通过以下步骤获取所述问题的解：

1)将子问题A_i划分为相应的对角子矩阵D_i及非对角子矩阵L_ij，并将右端项向量b进行划分，得到子右端项向量b_i；

2)为各对角子矩阵D_i及非对角子矩阵L_ij分配线程；

3)对于对角子矩阵，线程T_i基于对角子矩阵D_i进行求解；对于非对角子矩阵，利用公式y_i←D_i×(b_i-∑_0≤j＜iL_ij×y_j)进行求解；

4)拼接相应的子解，得到所述问题的解y。

进一步地，对于非对角子矩阵，通过以下步骤进行求解：

1)使用循环展开与SIMD向量化指令，并行执行普通矩阵-向量乘计算L_ij*y_j；

2)将各计算结果归约至线程T₀；

3)0号线程根据归约结果、对角子矩阵D_i及右端项向量段b_i进行回代求解，得到子解y_i。

进一步地，计算L_i(i-1)*y_(i-1)前，相应线程与线程T₀进行同步。

进一步地，线程之间的通信方法包括：RMA点对点通信。

一种存储介质，所述存储介质中存储有计算机程序，其中，所述计算机程序被设置为运行时执行以上所述方法。

一种电子装置，包括存储器和处理器，其中存储器存储执行以上所述方法的程序。

本发明取得了以下技术效果：

本发明实现了BLAS 1、2级函数的并行化。本发明设计了一种线程归约机制和一种线程通讯机制，解决了线程间的数据依赖问题。本发明还使用了循环变换和向量化技术优化计算。另外，本发明设计了一种自适应调优机制，根据输入问题的规模设置合适的线程数量，进一步提升了函数的性能。本发明的高性能BLAS 1、2级函数库与单核开源BLAS数学库GotoBLAS相比，平均加速比为22.37，最高加速比为65.47。

附图说明

图1为本发明面向SW26010-Pro处理器的1、2级BLAS函数库的高性能实现方法整体流程示意图；

图2为向量分段与核间数据映射示意图；

图3为普通矩阵分块与核间数据映射示意图；

图4为对称矩阵分块与核间数据映射示意图；

图5为三角矩阵分块与核间数据映射示意图；

图6为线程归约机制的示意图,其中(a)为线程行归约示意图，(b)为线程列归约示意图；

图7为TRSV的任务分块示意图；

图8为线程通讯机制的示意图；

图9为AXPY的任务分段示意图；

图10为GEMV的任务分块示意图；

图11为SYMV的任务分块示意图；

图12为本发明与开源库GotoBLAS的性能加速比。

具体实施方式

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

本发明的高性能实现方法，其特征包括：

特征一、根据输入问题的规模对矩阵或向量进行任务划分，产生若干子任务，并将各个子任务分配给每个线程。

特征二、提出一种基于RMA通信的线程归约机制和一种基于点对点同步的线程通讯机制。

特征三、使用循环变换和向量化技术优化计算。

特征四、提出一种自适应调优机制，为矩阵或向量的各个规模设置合适的线程数量。

进一步地，特征一包括：

如图2，对于向量，将其平均划分成若干向量段，将每个向量段依次映射给各个线程，图中的T₀，T₁，T₂，...T₆₃表示0号线程，1号线程，2号线程，…，63号线程；

如图3，对于普通矩阵，将其划分成若干小矩阵，将每个行块依次映射给各个线程，图中的T₀，T₁，T₂，...T₆₃表示0号线程，1号线程，2号线程，，...，63号线程；

如图4，对于对称矩阵，将其划分成若干小矩阵，将每个列块依次映射给各个线程，图中的T₀，T₁，T₂，...T₆₃表示0号线程，1号线程，2号线程，…，63号线程；

如图5，对于三角矩阵，将其划分成若干小矩阵，将对角块映射给0号线程，每个列块(对角块除外)映射给其它线程，图中的T₀，T₁，T₂，...T₆₃表示0号线程，1号线程，2号线程，...，63号线程。

进一步地，特征二包括：

如图6，从0号线程开始的任意个连续线程，采用RMA点对点通信进行归约，首先，位于同一行的线程进行行归约，归约目标为核组的第一列线程，而后，核组的第一列线程进行列归约，归约目标为0号线程，图中的T₀，T₁，T₂，...T₆₃表示0号线程，1号线程，2号线程，...，63号线程；

0号线程在完成当前操作后，会对1～63号中的某个线程发起一次点对点同步请求，与此同时，对应的线程在进行相应的操作前会响应这个同步请求。

进一步地，特征三包括：

对BLAS 1、2级函数的计算部分，使用了循环展开和SIMD向量化指令进行了优化。

进一步地，特征四包括：

根据向量规模将BLAS 1级函数分为四种类型，包括小规模、中等规模、大规模、超大规模。对于以上四种类型，分别启动8、16、32、64个线程。这里的小规模向量范围优选为[1024，4096]，中等规模向量范围优选为(4096，32768]，大规模向量范围优选为(32768，262144]，超大规模向量范围优选为(262144，+∞)；

根据矩阵规模将BLAS 2级函数分为两种类型，包括小规模、大规模。对于以上两种类型，分别启动16、64个线程。这里的小规模矩阵范围优选为[128*128，2048*2048]，大规模矩阵范围优选为(2048，+∞)。

下面以涉及下三角系数矩阵的线程方程组求解(TRSV)为例，它主要解决如下等式：A*x＝b，其中，A表示下三角矩阵，x表示待解未知向量，b表示右端项向量。具体实施步骤包括：

步骤一：根据矩阵A的规模，确定函数需要启动的线程数量。

步骤二：如图7所示，本发明按行对矩阵A、向量x和向量b进行任务划分，将每个行块视作一个子问题，共产生k个子问题。本发明对子问题进一步划分，每个子问题将产生对角子矩阵D_i(0≤i≤k-1)以及若干非对角子矩阵L_ij((0≤j＜i))，它们对应待解未知向量段x_i的解y_i和右端项向量段b_i。每个子问题完成如下运算：y_i←D_i×(b_i-∑_0≤j＜iL_ij×y_j)。

步骤三：本发明依次遍历每个子问题，将子问题对角部分的运算分给0号线程，非对角部分的运算按线程号依次分给其它线程。假设当前处理的是子问题i(0＜i≤k-1)，负责非对角部分的线程，使用循环展开和SIMD向量化指令并行执行普通矩阵-向量乘计算：L_ij*y_j，(0≤j＜i)，将计算结果归约至0号线程，0号线程根据归约结果、对角子矩阵D_i及右端项向量段b_i进行回代求解(back substitution)，得到y_i，并将y_i写回主存。例如当前处理的是子问题3，负责非对角线部分的线程1，线程2，线程3，并行执行普通矩阵-向量乘计算：L₃₀*x₀，L₃₁*x₁，L₃₂*x₂。将计算结果归约至0号线程，0号线程根据归约结果进行回代求解(back substitution)，得到x₃，并将x₃写回主存。

如图8所示，负责子矩阵L_i(i-1)的线程在计算前需要和0号线程同步，等待0号线程将y_(i-1)写回主存。L_ij的规模是128×128，L_ij*y_j采用两层循环实现，本发明将外层循环展开8次，增加了单次循环中乘加运算操作的数量，并在计算过程中使用SW26010-pro众核处理器硬件提供的浮点向量乘加指令加速了乘加运算。

步骤四：输出向量x的解y。

下面以标量向量乘(AXPY)为例，它的计算形式为：y＝α*x+y，其中，x和y表示向量，α表示标量。具体实施步骤包括：

步骤一：根据向量的规模确定函数需要启动的线程数量。

步骤二：如图9所示，本发明对向量x、待解未知向量y进行任务划分，将每个向量段视作一个子问题，共产生k个子问题。每个子问题完成如下运算：y_i←α×x_i+y_i。

步骤三：本发明依次遍历每个子问题，将子问题i分给i号线程。假设当前处理的是子问题i(0≤i≤63)，i号线程执行计算：α×x_i+y_i，得到y_i，并将y_i写回主存。

步骤四：输出向量y。

下面以普通矩阵向量乘(GEMV)为例，它的计算形式为：y＝α*A*x+β*y，其中，A表示普通矩阵，x和y表示向量，α和β表示标量。具体实施步骤包括：

步骤一：根据矩阵A的规模确定函数需要启动的线程数量。

步骤二：如图10所示，本发明按行对矩阵A、待解未知向量y进行任务划分，将每个行块视作一个子问题，共产生k个子问题。每个子问题完成如下运算：y_i←α×A_i×x+β×y_i。

步骤三：本发明依次遍历每个子问题，将子问题i分给i号线程。假设当前处理的是子问题i(0≤i≤63)，i号线程执行计算：α×A_i×x+β×y_i，得到y_i，并将y_i写回主存。

步骤四：输出向量y。

下面以涉及下三角矩阵的对称矩阵向量乘(SYMV)为例，它的计算形式为：y＝α*A*x+β*y，其中，A表示下三角对称矩阵，x和y表示向量，α和β表示标量。具体实施步骤包括：

步骤一：根据矩阵A的规模确定函数需要启动的线程数量。

步骤二：如图11所示，本发明按列对矩阵A进行任务划分，将每个列块视作一个子问题，共产生k个子问题。本发明对子问题进一步划分，每个子问题将产生对角子矩阵D_j(0≤j≤k-1)以及若干下三角子矩阵L_ij(i≥j)。每个子问题完成如下操作：对于对角子矩阵，将D_j的上三角部分用下三角的元素补齐，并计算：y_j←D_j×x_j+y_j；对于下三角子矩阵，计算：y_i←L_ij×x_j+y_i；对于下三角子矩阵的对称部分计算：y_j←L_ij×x_i+y_j。

步骤三：本发明依次遍历每个子问题，将子问题j分给j号线程。假设当前处理的是子问题j(0≤j≤(k-1))，j号线程执行操作：将D_j的上三角部分用下三角的元素补齐，并计算：D_j×x_j+y_j，得到y_j，并将y_j写回主存；计算：L_ij×x_j+y_i，得到y_i，并将y_i写回主存；计算：L_ij×x_i+y_j，得到y_j，并将y_j写回主存。

步骤四：输出向量y。

本实施例采用GotoBLAS数学库对本发明的性能加速效果进行验证。实施例选取的问题规模保证了两个版本的函数性能均达到了各自的最优值，选取的精度为实数双精度。图12为本发明与开源库GotoBLAS的性能加速比，从中可以看出本发明相对GotoBLAS的平均加速比为22.37，最高加速比为65.47。

本实施例仅针对SW26010-Pro众核处理器，将本发明的内容经简单变形后移植到其他平台，或者在没有对本发明的任务划分以及线程归约机制等进行创造性改进，抑或在执行计算阶段于本发明的基础上进行简单优化，本质上仍未脱离本发明所涵盖的内容，仍属于本发明保护的范畴。

本发明未详细阐述的部分属于本领域技术人员的公知技术。

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

Claims (3)

1.一种面向SW26010-Pro处理器的1、2级BLAS函数库的高性能实现方法，其步骤包括：

1)对问题进行任务划分，产生若干子问题，其中所述问题的结构包括向量、普通矩阵、对称矩阵或三角矩阵；其中，

在所述问题的结构为向量的情况下，将每一向量段视作一子问题x_i，其中0≤i≤k-1，k为子问题的数量；

在所述问题的结构为普通矩阵的情况下，将每个行块视作一子问题A_i；

在所述问题的结构为对称矩阵的情况下，每一列块视作一子问题A_j，其中j+1是矩阵的列编号，0≤j≤k-1；

在所述问题的结构为三角矩阵的情况下，每一行块视作一子问题A_i；

2)，将各子问题的运算分配给相应线程，以得到该线程的运算结果；其中，

在所述问题的结构是向量的情况下，所述将各子问题的运算分配给相应线程，以得到该线程的运算结果，包括：

将子问题x_i分配给相应线程T_i，以使线程T_i执行计算，得到线程T_i对应的解y_i＝α×x_i+y_i，α表示标量；

在所述问题的结构是普通矩阵的情况下，所述将各子问题的运算分配给相应线程，以得到该线程的运算结果，包括：

将子问题A_i分配给线程T_i，以使线程T_i执行计算，得到线程T_i对应的解y_i＝α×A_i×x+β×y_i，β表示标量，x表示向量；

在所述问题的结构是三角矩阵的情况下，所述将各子问题的运算分配给相应线程，以得到该线程的运算结果，包括：

对每一子问题A_i进行划分，获取对角子矩阵D_i及非对角子矩阵L_ij，0≤j＜i；

将对角子矩阵D_i的运算分给0号线程，将非对角子矩阵L_ij的运算按线程号i依次分给其它线程；

该其它线程使用循环展开和SIMD向量化指令并行执行普通矩阵-向量乘计算：L_ij*y_j，并将计算结果归约至0号线程；

0号线程根据归约结果、对角子矩阵D_i及右端项向量段b_i进行回代求解，得到0号线程对应的解y_i；

在所述问题的结构是对称矩阵的情况下，所述将各子问题的运算分配给相应线程，以得到该线程的运算结果，包括：

对每一子问题A_j进行划分，获取对角子矩阵D_j及非对角子矩阵L_ij，i≥j；

将子问题A_j分配给相应线程T_j；

应线程T_j将D_j的上三角部分用下三角的元素补齐，并计算：y_j←D_j×x_j+y_j；对于下三角子矩阵，计算：y_i←L_ij×x_j+y_i；对于下三角子矩阵的对称部分计算：y_j←L_ij×x_i+y_j；

3)拼接各线程的运算结果，获取所述问题的解。

2.一种存储介质，所述存储介质中存储有计算机程序，其中，所述计算机程序被设置为运行时执行权利要求1所述方法。

3.一种电子装置，包括存储器和处理器，所述存储器中存储有计算机程序，所述处理器被设置为运行所述计算机程序以执行如权利要求1所述方法。

Images