CN111859277B - 一种稀疏矩阵向量乘法向量化实现方法 - Google Patents

Abstract

本发明公开一种稀疏矩阵向量乘法向量化实现方法，步骤包括：步骤1.在DDR中分别构建第一矩阵、第二矩阵，按行读取待计算稀疏矩阵中非0元素数据值并按列存储在第一矩阵中，对应的列索引值按列存储在第二矩阵中；步骤2.在AM中配置数据缓冲区；步骤3.将待计算稠密向量、第二矩阵中数据分别从DDR传输到GSM；步骤4.从DDR中依次传输第一矩阵中数据到数据缓冲区；步骤5.读取第二矩阵中的列索引值，按照读取的列索引值读取待计算稠密向量的数据值，并传输到数据缓冲区；步骤6.对数据缓冲区中的数据执行向量化计算。本发明具有实现操作简单、资源利用率以及计算效率高、硬件开销小等优点。

Description

一种稀疏矩阵向量乘法向量化实现方法

技术领域

本发明涉及多核向量处理器技术领域，尤其涉及一种面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法。

背景技术

HPL(High Performance Linpack)是评测高性能计算机***性能的最重要基准测试程序，也是主要的评测指标。HPL求解的是大规模稠密线性方程组，具有良好的数据局部性，容易开发并行性和局部性。然而，随着高性能应用的不断发展，大量的科学计算应用，如大规模海洋天气预报，其应用特征与HPL并不相符。由此，自2013年以来，HPCG(HighPerformance Conjugate Gradient)被作为一种新的基准测试程序，用以补充高性能计算机***的性能评测。

HPCG作为主要的高性能计算机***性能测试基准之一，是采用区域分解(DomainDecomposition)、多重网格(Multigrid)和高斯-赛德尔预条件共轭梯度法(Preconditioned Conjugate Gradient)求解大型稀疏方程组，其特点是计算访存比低、内存访问不规则和细粒度的迭代计算。HPCG主要的计算时间耗费在ComputeSPMV和ComputeSYMGS两个函数上，其中ComputeSPMV即是稀疏矩阵向量乘法(Sparse MatrixVector Multiplication，SPMV)，而ComputeSYMGS与ComputeSPMV非常类似，主要也是SPMV算法。因此，SPMV是影响HPCG性能最主要的核心算法，HPCG性能的优化也即主要是SPMV算法的优化。

随着功耗和散热问题的日益突出，能耗逐步成为影响高性能计算***的越来越重要的因素，使得处理器的体系结构朝着多核方向发展，目前多核处理器已经占据了高性能计算市场的绝大部分。多核向量处理器是一种新颖的体系结构，在保持较低功耗的同时，具有强大的计算能力，适合加速大规模科学与工程计算。如图1所示，多核向量处理器通常包括多个核心，多核间通过片内的全局共享存储(Global Shared Memory，GSM)和片外的DDR存储***(DDR Memory System)实现数据共享和交换。每个向量处理器核心通常包括标量处理部件(Scalar Processing Unit，SPU)和向量处理部件(Vector Processing Unit，VPU)，SPU负责标量任务计算和流控，VPU负责向量计算。SPU包括可配置为Cache或SRAM模式的标量存储器SM。VPU包括VPE0，VPE1，......，VPEp多个向量计算单元，VPU支持向量数据的Load/Store，提供大容量的向量阵列存储器AM。

多核向量处理器计算性能高，可利于实现高效的SPMV算法优化，但是多核向量处理器的体系结构复杂，因此针对多核向量处理器的体系结构特点，研究高效的SPMV算法优化，对评测多核向量处理器的计算效率，发挥多核向量处理器的计算优势和提高应用程序的计算效率均具有重要的参考价值。

传统针对SPMV计算均是采用ELL格式存储数据，即按行存储稀疏矩阵的非0元素值，按行存储稀疏矩阵非0元素值的列索引值，按照列索引值取得maxnonzeros(每行中非零元素的最大个数)长度的向量数据值，该向量再与稀疏矩阵的非0元素值分别进行乘累加得到一个新的结果值。但是采用该类方式，一方面，maxnonzeros不一定是VPE数量的整数倍(例如maxnonzeros＝27，VPE数量＝16)，当maxnonzeros并不是VPE数量的整数倍时，会导致计算资源浪费；另一方面，向量的乘累加是在各个VPE之间进行，需要向量处理器支持混洗或归约计算，需要的硬件开销(面积、功耗等)大。因此亟需提供一种面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法，以使得能够充分发挥向量处理器体系结构的优势，提高向量处理器的资源利用率以及计算效率，同时减少硬件开销。

发明内容

本发明要解决的技术问题就在于：针对现有技术存在的技术问题，本发明提供一种实现操作简单、资源利用率以及计算效率高、硬件开销小的稀疏矩阵向量乘法向量化实现方法，能够充分发挥向量处理器体系结构的优势实现SPMV优化。

为解决上述技术问题，本发明提出的技术方案为：

一种面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法，步骤包括：

步骤1.按照TELL(Transposed ELLPACK)数据格式进行数据存储：在DDR中分别构建第一矩阵AV、第二矩阵AC，按行读取待计算稀疏矩阵A中非0元素数据值并按列存储在所述第一矩阵AV中，将各所述非0元素数据值所对应的列索引值按列存储在所述第二矩阵AC中；

步骤2.在AM中配置第一数据缓冲区BUFFER_AV1、第二数据缓冲区BUFFER_AX1；

步骤3.将待计算稠密向量x、所述第二矩阵AC中数据分别从DDR传输到GSM；

步骤4.从DDR中依次传输所述第一矩阵AV中数据到所述第一数据缓冲区BUFFER_AV1；

步骤5.读取所述第二矩阵AC中的列索引值，按照读取的所述列索引值从GSM中读取待计算稠密向量x的数据值，并传输到所述第二数据缓冲区BUFFER_AX1；

步骤6.向量处理器对所述第一数据缓冲区BUFFER_AV1、第二缓冲区BUFFER_AX1中的数据执行向量化计算，直至完成稀疏矩阵向量乘法Ax＝b中结果向量b的所有计算。

进一步的，所述步骤1中按照TELL数据格式进行数据存储时，还包括对所述第一矩阵AV中非0元素的个数小于maxnonzeros的列进行补0操作，并将所述第二矩阵AC中对应列的剩余元素的列索引值设置为指定负值步骤，其中maxnonzeros为所述待计算稀疏矩阵A的每行中非0元素的最大个数。

进一步的，所述步骤1的具体步骤包括：遍历所述待计算稀疏矩阵A的每一行；遍历过程中，取第i行为当前行，i＝0,1,2……,n，n为所述待计算稀疏矩阵A的行数，遍历所述待计算稀疏矩阵A的第i行，将第i行的每一个非0元素值依次存储在所述第一矩阵AV的第i列，将第i行的每一个非0元素值对应的列索引值依次存储在所述第二矩阵AC的第i列，其中若所述第二矩阵AV的第i列中非0元素个数小于maxnonzeros，则将该第i列的剩余元素的数据值设置为0，将所述第二矩阵AC中对应列的剩余元素的列索引值设为指定负值。

进一步的，所述第一数据缓冲区BUFFER_AV1的规模为maxnonzeros*(p*c*k)，所述第二数据缓冲区BUFFER_AX1的规模为maxnonzeros*(p*c*k)，其中p为向量处理器的VPE个数，c为VPE支持的FMAC单元数，AM的容量为s字节，计算的数据宽度为d字节，k为满足下式的最大整数：

(maxnonzeros*(p*c*k)+maxnonzeros*(p*c*k)+1*(p*c*k))*d<＝s/2

其中，maxnonzeros为所述待计算稀疏矩阵A的每行中非0元素的最大个数。

进一步的，所述步骤2中还包括配置第三数据缓冲区BUFFER_B1，所述第三数据缓冲区BUFFER_B1的规模为1*(p*c*k)，所述步骤6中每次计算得到的结果存储在所述第三数据缓冲区BUFFER_B1中，完成数据缓冲区的所有计算后，还包括将所述第三数据缓冲区BUFFER_B1的数据传输到DDR的步骤7。

进一步的，所述步骤2中具体采用AM双缓冲配置方式，即在AM的一半存储区域配置一份计算所需的所有数据缓冲区、在另一半存储区域再对应各数据缓冲区重新配置一份，以使得计算与数据传输重叠。

进一步的，所述步骤5中所述第一数据缓冲区BUFFER_AV1中数据值具体按入下规则计算得到：

若AC[i][j]>＝0，则BUFFER_AX1[i][j]＝x[AC[i][j]]；

若AC[i][j]<0，则BUFFER_AX1[i][j]＝0；

其中，i,j分别为行、列值。

进一步的，所述步骤6的具体步骤包括：

步骤6.1.初始化向量寄存器VRk的值为0；

步骤6.2.向量处理器分别从所述第一数据缓冲区BUFFER_AV1、所述第二数据缓冲区BUFFER_AX1的每一行加载数据到向量寄存器VRi、VRj，并按照VRk+＝VRi*VRj执行一次累加计算，直至完成所有行数据计算。

进一步的，该方法按照多核并行计算，所述步骤3中由多核向量处理器的每个核分别将所述第二矩阵AC中数据从DDR传输到GSM，所述步骤4中由多核处理器的每个核分别从DDR传输所述第一矩阵AV中数据到所述第一数据缓冲区BUFFER_AV1。

进一步的，所述步骤3中多核向量处理器的每个核具体根据各核id将所述第二矩阵AC中数据从DDR的地址address_AC传输到GSM中，各核负责计算任务的数据块地址address_AC为：

address_AC＝ADDRESS_AC+(num*LEN*q+id*LEN)*d

所述步骤4中多核向量处理器的每个核具体根据当前核id从DDR的地址address_AV传输所述第一矩阵AV中数据到数据缓冲区BUFFER_AV1，各核负责计算任务的数据块地址address_AV为：

address_AV＝ADDRESS_AV+(num*LEN*q+id*LEN)*d

其中，q为多核向量处理器的核数，ADDRESS_x为稠密向量x在DDR的起始地址，LEN＝p*c*k*2，p为向量处理器的VPE个数，c为VPE支持的FMAC单元数，AM的容量为s字节，计算的数据宽度为d字节，k为满足下式的最大整数：

(maxnonzeros*(p*c*k)+maxnonzeros*(p*c*k)+1*(p*c*k))*d<＝s/2

其中，maxnonzeros为所述待计算稀疏矩阵A的每行中非0元素的最大个数。

与现有技术相比，本发明的优点在于：

(1)本发明针对向量处理器的结构特性，构建一种存储稀疏矩阵数据的压缩数据存储格式TELL，向量处理器要进行SPMV计算时，按照该TELL存储格式按列存储稀疏矩阵的非0元素值以及非0元素值的列索引值，使得对任意的非0元素个数maxnonzeros，所有乘累加计算都是在同一个VPE上进行，不论待计算稀疏矩阵中maxnonzeros的大小是否为VPE数量的整数倍，各个VPE之间的数据均不需要进行混洗或者归约计算，因而无需向量处理器支持归约求和，可以有效减少硬件开销，同时大大提高向量处理器的资源利用率，能够充分发挥向量处理器体系结构的优势高效实现SPMV向量化。

(2)本发明按照TELL存储格式存储稀疏矩阵数据时，通过在数据存储中采用补0值的方式固定后续向量化计算的乘累加次数，基于固定的循环次数，可以便于采用软件流水高效实现，从而易于提升计算效率。

(3)本发明通过将存储有稀疏矩阵列索引值的第二矩阵和稠密向量从DDR传输到共享的且更快速的GSM中，能够尽可能多得让处理器从GSM中读取数据，减少从DDR读取数据的次数，从而减少整体计算时间。

附图说明

图1是典型多核向量处理器的结构示意图。

图2是本发明在具体应用实施例中构建的存储稀疏矩阵数据的压缩数据存储格式TELL的实现原理示意图。

图3是本实施例面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法的实现流程示意图。

图4是本发明在具体应用实施例中向量处理器进行SPMV向量化的具体流程示意图。

具体实施方式

以下结合说明书附图和具体优选的实施例对本发明作进一步描述，但并不因此而限制本发明的保护范围。

本实施例假定SPMV算法求解的问题为：Ax＝b，其中A为待计算稀疏矩阵，矩阵的行数为n，每行中非0元素的最大个数为maxnonzeros，x为待计算稠密向量，元素个数为n，b为结果向量，元素个数为n。如图3所示，本实施例面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法的步骤为：

步骤1.按照TELL(Transposed ELLPACK)数据格式进行数据存储：在DDR中分别构建第一矩阵AV、第二矩阵AC，按行读取待计算稀疏矩阵A中非0元素数据值并按列存储在第一矩阵AV中，将各非0元素数据值所对应的列索引值按列存储在第二矩阵AC中。

首先根据待计算稀疏矩阵A的行数和每行中非0元素的最大个数maxnonzeros，在DDR中分别构建两个连续存储的二维矩阵：第一矩阵AV和第二矩阵AC，两个矩阵的规模均为maxnonzeros*n，即行数为maxnonzeros、列数为n，其中，第一矩阵AV用于存储待计算稀疏矩阵A的元素的数据值，第二矩阵AC用于存储待计算稀疏矩阵A的元素列索引值。通过按照上述TELL数据格式进行数据存储，对任意的非0个数maxnonzeros和任意的VPE个数p，向量处理器均能高效的实现SPMV向量化。

进一步的，上述步骤S1按照TELL数据格式进行数据存储时，还包括对第一矩阵AV中非0元素的个数小于maxnonzeros的列进行补0操作，并将第二矩阵AC中对应列的剩余元素的列索引值设置为指定负值步骤，其中maxnonzeros为待计算稀疏矩阵A的每行中非0元素的最大个数。对于任意的稀疏矩阵A，每行的非0个数可能不一样，通过采用上述补0操作处理，使得在内核计算时能够固定乘累加次数(具体固定累加次数为maxnonzeros次)，可以便于采用软件流水高效实现，从而提高计算效率。

上述负值可以根据实际需求设置，本实施例具体配置为-1，也可以设置为任意其他负值。

在具体应用实施例中，上述步骤1的具体步骤包括：

遍历待计算稀疏矩阵A的每一行；

遍历过程中，取第i行为当前行，i＝0,1,2……,n，n为待计算稀疏矩阵A的行数，遍历待计算稀疏矩阵A的第i行，将第i行的每一个非0元素值依次存储在第一矩阵AV的第i列，将第i行的每一个非0元素值对应的列索引值依次存储在第二矩阵AC的第i列，其中若第二矩阵AV的第i列中非0元素个数小于maxnonzeros，则将该第i列的剩余元素的数据值设置为0，将第二矩阵AC中对应列的剩余元素的列索引值设为指定负值。

上述步骤完毕以后，得到存储有稀疏矩阵A的元素的数据值的第一矩阵AV，以及存储有对应的列索引值的第二矩阵AC，两个矩阵的规模均为maxnonzeros*n。

步骤2.在AM中配置第一数据缓冲区BUFFER_AV1、第二数据缓冲区BUFFER_AX1。

上述第一数据缓冲区BUFFER_AV1用于缓存第一矩阵AV中的数据，第二数据缓冲区BUFFER_AX1用于缓存依据第二矩阵AC中列索引值从稠密矩阵x中读取的数据值。

设向量处理器的VPE个数为p，VPE支持的FMAC单元数为c，AM的容量为s字节，计算的数据宽度为d字节。上述数据缓冲区BUFFER_AV1、BUFFER_AX1的规模具体分别配置为：maxnonzeros*(p*c*k)，maxnonzeros*(p*c*k)，其中k为满足下式的最大整数：

(maxnonzeros*(p*c*k)+maxnonzeros*(p*c*k)+1*(p*c*k))*d<＝s/2 (1)

上述步骤2中还包括配置第三数据缓冲区BUFFER_B1，以用于缓存后续向量化计算时的结果，该第三数据缓冲区BUFFER_B1的规模具体配置为1*(p*c*k)。

步骤3.将待计算稠密向量x、第二矩阵AC中数据分别从DDR传输到GSM。

将与上述数据缓冲区规模相对应的第二矩阵AC、稠密向量x分别从DDR传输到GSM中。通过将存储有稀疏矩阵列索引值的第二矩阵AC和稠密向量x从DDR传输到共享的且更快速的GSM中，能够尽可能多得让处理器从GSM中读取数据，减少从DDR读取数据的次数，从而大大减少整体计算时间。

步骤4.从DDR中依次传输第一矩阵AV中数据到第一数据缓冲区BUFFER_AV1。

具体按照数据缓冲区的规模，从DDR中依次传输第一矩阵AV中的数据到第一数据缓冲区BUFFER_AV1，以将待计算稀疏矩阵中非0元素的数据值传输至数据缓冲区等待计算。

步骤5.读取第二矩阵AC中的列索引值，按照读取的列索引值从GSM中读取待计算稠密向量x的数据值，并传输到第二数据缓冲区BUFFER_AX1。

具体按照数据缓冲区的规模，依次从GSM读取步骤4中第一AV对应的AC列索引值，按照读取的AC列索引值数据从GSM读取稠密向量x的数据值，传输到第二数据缓冲区BUFFER_AX1，以按照AC索引将稠密向量x中数据传输到数据缓冲区等待计算。

上述具体按照如下规则读取待计算稠密向量x的数据值：

若AC[i][j]>＝0，则BUFFER_AX1[i][j]＝x[AC[i][j]]；即AC[i][j]的值作为x向量的索引值，得到该索引值的x向量的元素值，例如，AC[i][j]＝5，则x[AC[i][j]]表示x向量的第5个元素值(0开始)，即x[5]。

若AC[i][j]<0，则BUFFER_AX1[i][j]＝0；

其中，i,j分别为行、列值。

步骤6.向量处理器对第一数据缓冲区BUFFER_AV1、第二缓冲区BUFFER_AX1中的数据执行向量化计算，直至完成稀疏矩阵向量乘法Ax＝b中结果向量b的所有计算。

上述执行向量化计算的具体步骤包括：

S61.初始化向量寄存器VRk的值为0；

S62.向量处理器分别从第一数据缓冲区BUFFER_AV1、第二数据缓冲区BUFFER_AX1的每一行加载数据到向量寄存器VRi、VRj，并按照VRk+＝VRi*VRj执行一次累加计算，直至完成所有行数据计算。

上述每次计算得到的结果存储在第三数据缓冲区BUFFER_B1中，具体固定执行maxnonzeros次累加计算后，将向量寄存器VRk的值存入BUFFER_B1。

步骤7.向量处理器完成数据缓冲区的所有计算后，将第三数据缓冲区BUFFER_B1的数据传输到DDR。

重复循环执行上述步骤2-7，直到完成向量b的所有计算。循环执行过程中，稠密向量x的起始位置为上一轮循环第二矩阵AC中的最大列索引值的下一个索引值，第一次是从0位置开始。

上述步骤2中设置数据缓冲区只用到了AM容量的一半，本实施例进一步采用AM双缓冲配置方式，即在AM的一半存储区域配置一份计算所需的所有数据缓冲区、在另一半存储区域再对应各数据缓冲区重新配置一份，以使得计算与数据传输重叠。具体可在AM中另一半存储区域中配置与上述数据缓冲区BUFFER_AV1，BUFFER_AX1和BUFFER_B1一一对应的3个数据缓冲区：BUFFER_AV2、BUFFER_AX2和BUFFER_B2，再对上述步骤2-7采用双缓冲策略计算，使得计算与数据传输重叠，可以进一步提高计算效率。

本实施例上述方法采用多核并行计算，以进一步提高计算效率。当采用多核并行计算时，步骤3由多核向量处理器的每个核分别将第二矩阵AC中数据从DDR传输到GSM，步骤S4中由多核处理器的每个核分别从DDR传输第一矩阵AV中数据到第一数据缓冲区BUFFER_AV1。

设多核向量处理器的核数为q，编号依次为0,1,2,...,q-1。设二维矩阵AV和AC在DDR的起始地址分别为ADDRESS_AV，ADDRESS_AC，稠密向量x在DDR的起始地址是ADDRESS_x。根据上述步骤2计算得到的k值，令LEN＝p*c*k*2。

设本次计算循环次数为num，则采用多核并行计算时：

步骤3中多核向量处理器的每个核根据各核id将与上述数据缓冲区规模相对应的存储稀疏矩阵列索引值的二维矩阵AC从DDR的地址address_AC传输到GSM中，将稠密向量x从DDR传输到GSM中；各核负责计算任务的数据块地址address_AC为：

address_AC＝ADDRESS_AC+(num*LEN*q+id*LEN)*d (2)

其中，稠密向量x的起始位置为上一轮循环各核中的AC的最大列索引值的下一个索引值，第一次是从0位置开始。

步骤4中按照数据缓冲区的规模，多核向量处理器的每个核根据各核id从DDR的地址address_AV传输AV数据到数据缓冲区BUFFER_AV1；各核负责计算任务的数据块地址address_AV为：

address_AV＝ADDRESS_AV+(num*LEN*q+id*LEN)*d (3)

本实施例针对向量处理器的结构特性，构建一种存储稀疏矩阵数据的压缩数据存储格式TELL，向量处理器要进行SPMV计算时，先按照该TELL存储格式按列存储稀疏矩阵的非0元素值，以及非0元素值的列索引值，使得对任意的非0元素个数maxnonzeros，乘累加都是在同一个VPE上进行，不论待计算稀疏矩阵中maxnonzeros的大小是否为VPE数量的整数倍，各个VPE之间的数据均不需要进行混洗或归约计算，无需向量处理器支持归约求和，可以有效减少硬件开销；且由于稀疏矩阵的行数n通常较大，可以设置为VPE数量的整数倍，而即便n不是VPE数量的整数倍，也仅有尾数处理这一极少部分没有利用全部的VPE计算资源，可以大大提高资源利用率。举例来说，n＝10007，VPE数量＝16，10007＝625*16+7，那么有625次循环计算是完全使用VPE的全部计算资源，最后7列的处理使用部分VPE计算资源，利用率可以达99.93％

对于任意的稀疏矩阵A，每行的非0个数可能不一样，传统是采用精确记录每行的非0个数k的方式以便在向量计算时精确进行k次乘累加计算，由于每行的非0个数k不一致，会导致向量计算时，每个VPE的计算长度可能不同，不方便软件流水，反而不利于提升计算效率。本实施例按照TELL存储格式存储稀疏矩阵A数据时，通过在数据存储中采用补0值的方式固定后续向量化计算的乘累加次数maxnonzeros，基于固定的循环次数，可以方便采用软件流水高效实现，易于提升计算效率。

在具体应用实施例中构建的存储稀疏矩阵A数据的压缩数据存储格式TELL如图2所示，其中稀疏矩阵A的行数和列数n＝8，每行中非0元素的最大个数maxnonzeros＝3。首先在DDR中分别构建两个连续存储的二维矩阵AV和AC，二维矩阵的规模均为3*8；遍历稀疏矩阵A的各行进行数据存储：

遍历稀疏矩阵A的第0行；得到非0元素数据1,3,7，依次存入AV的第0列，对应的列索引值1,3,5，依次存储在AC的第0列；

遍历稀疏矩阵A的第1行；得到非0元素数据2,8，依次存入AV的第1列，对应的列索引值0,4，依次存储在AC的第1列；非0元素个数2小于最大非0个数maxnonzeros＝3，AV第1列剩余元素值设为0，AC第1列剩余元素值设为-1；

遍历稀疏矩阵A的第2行；得到非0元素数据5,9，依次存入AV的第2列，对应的列索引值1,5，依次存储在AC的第2列；非0元素个数2小于最大非0个数maxnonzeros＝3，AV第2列剩余元素值设为0，AC第2列剩余元素值设为-1；

遍历稀疏矩阵A的第3行；得到非0元素数据6,7，依次存入AV的第3列，对应的列索引值2,6，依次存储在AC的第3列；非0元素个数2小于最大非0个数maxnonzeros＝3，AV第3列剩余元素值设为0，AC第3列剩余元素值设为-1；

继续上述过程；

遍历稀疏矩阵A的第6行；得到非0元素数据5,9，依次存入AV的第6列，对应的列索引值1,5，依次存储在AC的第6列；非0元素个数2小于最大非0个数maxnonzeros＝3，AV第6列剩余元素值设为0，AC第6列剩余元素值设为-1；

遍历稀疏矩阵A的第7行；得到非0元素数据6,7，依次存入AV的第7列，对应的列索引值2,6，依次存储在AC的第7列；非0元素个数2小于最大非0个数maxnonzeros＝3，AV第7列剩余元素值设为0，AC第7列剩余元素值设为-1；

上述计算完毕，得到存储稀疏矩阵A的元素的数据值的第一矩阵AV，以及存储对应的列索引值的第二矩阵AC，二维矩阵的规模均为3*8。

在具体应用实施例中向量处理器实现SPMV计算如图4所示，具体实施步骤为：

采用如图2中稀疏矩阵A，稠密向量x的数据值为{1,2,3,4,5,6,7,8}；

将第一矩阵AV从DDR传输到AM中得到第一数据缓冲区BUFFER_AV1；

将第二矩阵AC和稠密向量x分别从DDR传输到GSM中；

从GSM中读取矩阵AC的列索引值，根据列索引值从GSM中读取稠密向量x数据值，传输到AM中得到第二数据缓冲区BUFFER_AX1；

将第一数据缓冲区BUFFER_AV1和第二数据缓冲区BUFFER_AX1中每一行向量分别进行乘累加，3次乘累加得到结果向量b。

本发明上述面向向量处理器的稀疏矩阵向量乘法向量化实现方法，对任意的非0元素个数maxnonzeros和任意的VPE个数p，不需要向量处理器支持归约求和，向量处理器都能高效的实现SPMV向量化，在内核计算时能固定乘累加次数，可以减少硬件开销，同时可以方便的采用软件流水高效实现，同时采用共享数据方式能够显著减少从DDR读取数据的次数，大幅提高计算效率。

上述只是本发明的较佳实施例，并非对本发明作任何形式上的限制。虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明。因此，凡是未脱离本发明技术方案的内容，依据本发明技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均应落在本发明技术方案保护的范围内。

Claims (10)

1.一种面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法，其特征在于，步骤包括：

步骤1.按照TELL数据格式进行数据存储：在DDR中分别构建第一矩阵AV、第二矩阵AC，按行读取待计算稀疏矩阵A中非0元素数据值并按列存储在所述第一矩阵AV中，将各所述非0元素数据值所对应的列索引值按列存储在所述第二矩阵AC中；

步骤2.在AM中配置第一数据缓冲区BUFFER_AV1、第二数据缓冲区BUFFER_AX1；

步骤3.将待计算稠密向量x、所述第二矩阵AC中数据分别从DDR传输到GSM；

步骤4.从DDR中依次传输所述第一矩阵AV中数据到所述第一数据缓冲区BUFFER_AV1；

步骤5.读取所述第二矩阵AC中的列索引值，按照读取的所述列索引值从GSM中读取待计算稠密向量x的数据值，并传输到所述第二数据缓冲区BUFFER_AX1；

步骤6.向量处理器对所述第一数据缓冲区BUFFER_AV1、第二缓冲区BUFFER_AX1中的数据执行向量化计算，直至完成稀疏矩阵向量乘法Ax＝b中结果向量b的所有计算。

2.根据权利要求1所述的面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法，其特征在于：所述步骤1中按照TELL数据格式进行数据存储时，还包括对所述第一矩阵AV中非0元素的个数小于maxnonzeros的列进行补0操作，并将所述第二矩阵AC中对应列的剩余元素的列索引值设置为指定负值步骤，其中maxnonzeros为所述待计算稀疏矩阵A的每行中非0元素的最大个数。

3.根据权利要求2所述的面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法，其特征在于，所述步骤1的具体步骤包括：遍历所述待计算稀疏矩阵A的每一行；遍历过程中，取第i行为当前行，i＝0,1,2……,n，n为所述待计算稀疏矩阵A的行数，遍历所述待计算稀疏矩阵A的第i行，将第i行的每一个非0元素值依次存储在所述第一矩阵AV的第i列，将第i行的每一个非0元素值对应的列索引值依次存储在所述第二矩阵AC的第i列，其中若所述第二矩阵AV的第i列中非0元素个数小于maxnonzeros，则将该第i列的剩余元素的数据值设置为0，将所述第二矩阵AC中对应列的剩余元素的列索引值设为指定负值。

4.根据权利要求2所述的面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法，其特征在于：所述第一数据缓冲区BUFFER_AV1的规模为maxnonzeros*(p*c*k)，所述第二数据缓冲区BUFFER_AX1的规模为maxnonzeros*(p*c*k)，其中p为向量处理器的VPE个数，c为VPE支持的FMAC单元数，AM的容量为s字节，计算的数据宽度为d字节，k为满足下式的最大整数：

(maxnonzeros*(p*c*k)+maxnonzeros*(p*c*k)+1*(p*c*k))*d<＝s/2

其中，maxnonzeros为所述待计算稀疏矩阵A的每行中非0元素的最大个数。

5.根据权利要求4所述的面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法，其特征在于，所述步骤2中还包括配置第三数据缓冲区BUFFER_B1，所述第三数据缓冲区BUFFER_B1的规模为1*(p*c*k)，所述步骤6中每次计算得到的结果存储在所述第三数据缓冲区BUFFER_B1中，完成数据缓冲区的所有计算后，还包括将所述第三数据缓冲区BUFFER_B1的数据传输到DDR的步骤7。

6.根据权利要求1～5中任意一项所述的面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法，其特征在于：所述步骤2中具体采用AM双缓冲配置方式，即在AM的一半存储区域配置一份计算所需的所有数据缓冲区、在另一半存储区域再对应各数据缓冲区重新配置一份，以使得计算与数据传输重叠。

7.根据权利要求1～5中任意一项所述的面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法，其特征在于，所述步骤5中所述第一数据缓冲区BUFFER_AV1中数据值具体按入下规则计算得到：

若AC[i][j]>＝0，则BUFFER_AX1[i][j]＝x[AC[i][j]]；

若AC[i][j]<0，则BUFFER_AX1[i][j]＝0；

其中，i,j分别为行、列值。

8.根据权利要求1～5中任意一项所述的面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法，其特征在于，所述步骤6的具体步骤包括：

步骤6.1.初始化向量寄存器VRk的值为0；

步骤6.2.向量处理器分别从所述第一数据缓冲区BUFFER_AV1、所述第二数据缓冲区BUFFER_AX1的每一行加载数据到向量寄存器VRi、VRj，并按照VRk+＝VRi*VRj执行一次累加计算，直至完成所有行数据计算。

9.根据权利要求1～5中任意一项所述的面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法，其特征在于，该方法按照多核并行计算，所述步骤3中由多核向量处理器的每个核分别将所述第二矩阵AC中数据从DDR传输到GSM，所述步骤4中由多核处理器的每个核分别从DDR传输所述第一矩阵AV中数据到所述第一数据缓冲区BUFFER_AV1。

10.根据权利要求9所述的面向多核向量处理器的稀疏矩阵向量乘法向量化实现方法，其特征在于，所述步骤3中多核向量处理器的每个核具体根据各核id将所述第二矩阵AC中数据从DDR的地址address_AC传输到GSM中，各核负责计算任务的数据块地址address_AC为：

address_AC＝ADDRESS_AC+(num*LEN*q+id*LEN)*d

其中，稠密向量x的起始位置为上一轮循环各核中的AC的最大列索引值的下一个索引值，第一次是从0位置开始；

所述步骤4中多核向量处理器的每个核具体根据当前核id从DDR的地址address_AV传输所述第一矩阵AV中数据到数据缓冲区BUFFER_AV1，各核负责计算任务的数据块地址address_AV为：

address_AV＝ADDRESS_AV+(num*LEN*q+id*LEN)*d

其中，q为多核向量处理器的核数，LEN＝p*c*k*2，p为向量处理器的VPE个数，c为VPE支持的FMAC单元数，AM的容量为s字节，计算的数据宽度为d字节，k为满足下式的最大整数：

(maxnonzeros*(p*c*k)+maxnonzeros*(p*c*k)+1*(p*c*k))*d<＝s/2

其中，maxnonzeros为所述待计算稀疏矩阵A的每行中非0元素的最大个数。

Images