阅读说明：本技术 一种适用于gpu硬件的矩阵存储与计算方法 (Matrix storage and calculation method suitable for GPU hardware ) 是由 邵雪 王晓光 周振亚 于 2019-09-11 设计创作，主要内容包括：一种适用于GPU硬件的矩阵存储与计算方法,包括以下步骤：1)存储矩阵的行数、列数、非零元素、每个元素是否非零的标志,以及起始非零元素所在位置；2)通过GPU访问矩阵元素,取得矩阵元素是否非零以及非零元素的值,设置矩阵非零元素的值；3)使用GPU进行矩阵运算。本发明的适用于GPU硬件的矩阵存储与计算方法,能实现GPU硬件下多线程高速访问矩阵中的任意元素,从而大幅提升了GPU中矩阵计算的速度。(a matrix storage and calculation method suitable for GPU hardware comprises the following steps: 1) storing the row number, the column number, the nonzero elements, a mark of whether each element is nonzero or not and the position of the initial nonzero element; 2) accessing matrix elements through a GPU (graphics processing Unit), acquiring whether the matrix elements are nonzero and the values of the nonzero elements, and setting the values of the nonzero elements of the matrix; 3) and performing matrix operation by using the GPU. The matrix storage and calculation method suitable for GPU hardware can realize multi-thread high-speed access to any element in the matrix under the GPU hardware, thereby greatly improving the matrix calculation speed in the GPU.)

一种适用于GPU硬件的矩阵存储与计算方法

技术领域

本发明涉及GPU硬件高性能计算领域，具体涉及GPU硬件对矩阵乘法计算及LU分解的高性能计算技术领域，特别涉及一种适用于GPU硬件的矩阵存储与计算方法。

背景技术

近年来高性能计算中矩阵运算的规模越来越大，所需要的计算能力也越来越强，传统CPU构架受限于功耗瓶颈难以进一步提高性能，无法胜任计算的需求。相比之下GPU具有计算资源充足、数据访问带宽高的优势，理想情况下相比CPU可以加速十数倍。但矩阵分解存在高相关性导致算法优化难度很大，GPU应用进展较为缓慢。

发明内容

为了解决现有技术存在的不足，本发明的目的在于提供一种适用于GPU硬件的矩阵存储与计算方法，充分利用GPU硬件的特性，实现矩阵的高性能计算。

为实现上述目的，本发明提供的适用于GPU硬件的矩阵存储与计算方法，包括以下步骤：

1)存储矩阵的行数、列数、非零元素、每个元素是否非零的标志，以及起始非零元素所在位置；

2)通过GPU访问矩阵元素,取得矩阵元素是否非零以及非零元素的值，设置矩阵非零元素的值；

3)使用GPU进行矩阵运算。

进一步地，所述步骤1)进一步包括：

存储矩阵的行数与列数；

按照矩阵行或列的顺序依次存储矩阵的非零元素至第一数组；

存储每行或列中每个元素是否非零的标志至第二数组；

存储每行或列的起始非零元素在所述第一数组中的位置至第三数组。

进一步地，在所述按照矩阵行或列的顺序依次存储矩阵的非零元素至第一数组的步骤还包括，根据矩阵的非零元总数目确定第一数组的大小。

进一步地，在所述存储每行或列中每个元素是否非零的标志至第二数组的步骤还包括，根据矩阵的行或列数确定第二数组的大小，将行或列的矩阵元素的非零标志按照从低到高的顺序分别连续存储至第二数组中，其中每一个bit对应一个矩阵元素，bit值为1代表对应的矩阵元素非零。

进一步地，所述步骤2)进一步包括以下步骤：

根据待读取数据的位置信息，取得其非零标志位；

读取待读取数据所在行或列的首个非零元素的位置，记作第一位置；

计算待读取数据与其所在行或列的首个非零元素之间的位置差，记作第二位置；

计算待读取数据在所述第一数组中的位置，记作第三位置，第三位置＝第一位置+第二位置；

根据第三位置在第一数组中读取待读取数据的值。

进一步地，所述步骤3)包括矩阵加法、矩阵减法、矩阵乘法和矩阵LU分解算法。

进一步地，所述矩阵加法包括以下步骤：

判断参与加法运算的两个矩阵中的两个对应数据是否非零；

如果所述两个对应数据均为零元素，则两个对应数据的加运算或减运算结果为零；

如果所述两个对应数据中只有一个数据不为零，则两个对应数据的加运算或减运算的结果为非零数据的正值或负值；

如果所述两个对应数据均不为零，则两个数据的加运算或减运算为所述两个对应数据的加和或相减。

进一步地，所述矩阵乘法包括以下步骤：

首先对C_ij进行判断，如果C_ij是非零元素则结束，否则继续；

初始化变量v＝0；

遍历A矩阵第i行从1到C_A的元素a_ik与B矩阵第j列从1到R_B的元素b_kj，若a_ik与b_kj均为非零元素则v＝v+a_ik×b_kj；

得到结果C_ij＝v。

进一步地，所述矩阵的LU分解算法，包括以下步骤：

将压缩矩阵数据读取到稠密矩阵D中，其中，如果a_ij为非零元素d_ij＝a_ij，否则d_ij＝0，依次同步GPU的所有线程，保证D矩阵中数据一致性；

按照从1到R或C的顺序遍历分解D矩阵的第k行第k列，重复以下步骤：

取D矩阵第k行从k到C列的元素为U矩阵的第k行结果，其中，如果u_kj为非零元素则u_kj＝d_kj，j的范围为从k+1到C；

取D矩阵第k列从k+1到R列的元素除以d_kk得到L矩阵的第k列结果，其中，如果l_ik为非零元素则l_ik＝d_ik/d_kk，i的范围为从k+1到R；

更新D矩阵第k行k列右下侧的所有元素，其中，如果l_ik与u_kj均为非零元素则d_ij＝d_ij-l_ik×u_kj，i的范围为从k+1到R，j的范围为从k+1到C；

同步GPU的所有线程，保证D矩阵中所有数据的一致性。

为实现上述目的，本发明还提供一种计算机可读存储介质，其上存储有计算机指令，所述计算机指令运行时执行上述的适用于GPU硬件的矩阵存储与计算方法的步骤。

有益效果：本发明的适用于GPU硬件的矩阵存储与计算方法，使用非零标志计算矩阵非零元素在存储空间中的位置，能实现GPU硬件下多线程高速访问矩阵中的任意元素，提高了GPU硬件对矩阵元素的访问效率，从而大幅提升了GPU硬件上矩阵的计算速度。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，并与本发明的实施例一起，用于解释本发明，并不构成对本发明的限制。在附图中：

图1为根据本发明的适用于GPU硬件的矩阵存储与计算方法的流程图；

图2为根据本发明的实施方式的分解的矩阵存储数据的示意图；

图3为根据本发明的实施方式的计算方法分解的矩阵LU分解的流程图。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

图1为根据本发明的适用于GPU硬件的矩阵存储与计算方法的流程图，下面将参考图1，对本发明的适用于GPU硬件的矩阵存储与计算方法进行详细描述。

在步骤11，存储数据。

在该步骤中，存储数据进一步包括以下步骤：

(111)存储矩阵的行数与列数；

(112)按照矩阵行(列)的顺序依次存储矩阵的非零元素；

(113)存储每行(列)中每个元素是否非零的标志；

(114)存储每行(列)的起始非零元素所在位置。

图2为根据本发明的实施方式的分解的矩阵存储数据的示意图，以下结合图2对步骤(111)-步骤(114)进行说明。

首先，在步骤(111)，存储矩阵的行数与列数。

在图2所示的实施方式中，矩阵的行数为5，列数为6，对其进行存储。

在步骤(112)，按照矩阵行(列)的顺序依次存储矩阵的非零元素。

在该步骤中，对于已有的矩阵根据需求确定其数据的存储顺序，以下说明中以行优先方式为例进行说明。根据矩阵的非零元总数目S为其准备矩阵元的存储空间M[S](如图2中101所示)。然后，按照矩阵行(列)的顺序依次存储矩阵的非零元素。

具体结合图2，矩阵的非零元总数目S＝12，需准备的矩阵元的存储空间为M[12]。然后，将5×6矩阵中的非零元素依次存储至M[12]中。

在步骤(113)，存储每行(列)中每个元素是否非零的标志。

在该步骤中，根据矩阵的行(列)数R为其准备存储非零元标志的数组F[R](如图2中102所示)，这里F[i-1]是第i行(列)中所有矩阵元素的非零标志按照从低到高的顺序连续存储，每一个bit对应一个矩阵元素，bit值为1代表这个矩阵元素非零。

对于列(行)数不大于32的矩阵F[i]所需的存储空间为一个32位的int变量，对于列(行)数不大于64的矩阵F[i]所需的存储空间为一个64位的int变量，超过64列(行)的矩阵F[i]所需的存储空间为ceil(矩阵列数/64)个64位的int变量，ceil为向上取整运算，即当矩阵列(行)数大于64时，每一行(列)的状态连续存储在多个64位int型变量中。

在步骤(114)，存储每行(列)的起始非零元素所在位置。

在该步骤中，所指的位置为每行(列)中第一个非零元素在(112)所指的连续存储空间M[S]中的位置序数，在其最后额外添加一个数据空间存储整个矩阵的非零元素总数，也即是矩阵最后一行(列)最后一个非零元素在(112)中的M[S]中的位置序数加一。

具体地，在该步骤中，根据矩阵的行数R为其准备存储行(列)起始位置的数组P[R+1](如图2中103所示)，P[0]＝0即M的起始位置，P[i]＝P[i-1]+第i行(列)的非零元素个数，P[R]＝矩阵的非零元素个数。

因此，一个矩阵A按行进行数据存储时，可以表示为以下数据的集合：

·矩阵的行数R、列数C、非零元总数S

·非零元数据数组M，长度为S

·非零元标志数组F，长度为R×ceil(C/64)

·行首位置数组P，长度为R+1

若矩阵列数不大于32，非零元标志数组的长度为R。

在步骤12，读取数据。

在该步骤中，基于步骤11中的数据存储进行数据读取。具体地，该步骤进一步包括以下步骤：

(121)读取矩阵存储数据的数组F中第i个数据的第j个bit是否为true，判断矩阵第i行(列)j列(行)元素是否非零；

(122)读取矩阵存储数据的数组P中第i个数据取得第i行(列)第一个非零元素在矩阵存储数据的数组M中的位置P_i；

(123)计算矩阵存储数据的数组F中第i个数据的第j-1个bit至第1个bit之和取得第i行(列)j列(行)元素相对本行(列)第一个非零元素在矩阵存储数据的数组M的相对位置Q_i；

(124)P_i与Q_i相加即可得到第i行(列)j列(行)元素在矩阵存储数据的数组M中的位置，然后进行取值、赋值操作。

具体地，在步骤(121)，判断矩阵元素a_ij是否非零。

取出非零元标志数组F中的第(i-1)×ceil(C/64)+(j-1)/64个数据f，

计算j-1对64取模的结果r，

取数据f的第r位，为1则a_ij是非零元素，为0则a_ij是零；

对于不大于32列的矩阵则简化为：

取非零元标志数组F中的第(i-1)个数据f，

取数据f的第j-1位，为1则a_ij是非零元素，为0则a_ij是零；

在步骤(122)，读取矩阵第i行的第一个非零元素位置p为P[i-1]。

在步骤(123)，计算矩阵元素a_ij相对于第i行的第一个非零元素位置偏移。

初始化位置偏移q为0，

遍历第i行矩阵从1到j-1的元素a_ik，若a_ik为非零元素则q＝q+1；

对于j达到64的情况，可以使用GPU的指令一次性取出前64位的非零元标志数据bit并求和以提高速度。

在步骤(124)，读取矩阵元素a_ij＝M[p+q]，写入矩阵元素M[p+q]＝a_ij(仅在a_ij为非零元素时有效)，进行赋值操作。

在步骤13，进行数据计算。

在该步骤中，基于步骤11和步骤12的数据存储与数据读取方法可以实现矩阵的加(减)法、乘法与LU分解运算。

对于矩阵的加(减)法，假设参与运算的两个矩阵分别为矩阵A和矩阵B，分别判断A和B矩阵位于第i行j列处的元素a_ij和b_ij是否非零，结果分别为f_A和f_B：

(1)若f_A和f_B均为1，读取矩阵A和B位于第i行j列处的元素值进行加(减)运算后作为结果矩阵C第i行j列的元素值；

(2)若f_A为1而f_B为0，取矩阵A位于第i行j列处的元素值作为结果矩阵C第i行j列的元素值；

(3)若f_A为0而f_B为1，取矩阵B位于第i行j列处的正(负)元素值作为结果矩阵C第i行j列的元素值；

(4)若f_A和f_B均为0，结果矩阵C第i行j列的元素值也为零。

对于矩阵A与B的乘法，要求矩阵A的列数与矩阵B的行数相等，即C_A与R_B相等，可以通过如下的算法实现：

首先对C_ij进行判断，如果C_ij是非零元素则结束，否则继续；

初始化变量v＝0；

遍历A矩阵第i行从1到C_A的元素a_ik与B矩阵第j列从1到R_B的元素b_kj，若a_ik与b_kj均为非零元素则v＝v+a_ik×b_kj；

得到结果C_ij＝v。

矩阵的加减法与乘法运算每个位置的元素都是多线程安全的，可以利用GPU运算中的多线程并行算法，每个线程计算一个位置的结果。

矩阵的LU分解运算则需要使用共享内存与线程同步技术进行线程之间的数据访问。

图3为根据本发明的实施方式的计算方法分解的矩阵LU分解的流程图，下面将参考图3对本发明的矩阵LU分解方法进行详细描述：

(1)将本发明所描述的压缩矩阵(矩阵的行列数应一致)数据读取到一段完整且连续的稠密矩阵D中，该稠密矩阵存储空间应为共享内存以保证GPU任意线程可以访问，空间为R×C：

如果a_ij为非零元素d_ij＝a_ij，否则d_ij＝0；

同步GPU的所有线程，保证D矩阵中数据一致性；

(2)按照从1到R(C)的顺序遍历分解D矩阵的第k行第k列，重复下面的步骤(3)至(6)；

(3)取D矩阵第k行从k到C列的元素为U矩阵的第k行结果：

如果u_kj为非零元素则u_kj＝d_kj，j的范围为从k+1到C；

(4)取D矩阵第k列从k+1到R列的元素除以d_kk得到L矩阵的第k列结果：

如果l_ik为非零元素则l_ik＝d_ik/d_kk，i的范围为从k+1到R；

(5)更新D矩阵第k行k列右下侧的所有元素：

如果l_ik与u_kj均为非零元素则d_ij＝d_ij-l_ik×u_kj，i的范围为从k+1到R，j的范围为从k+1到C；

(6)同步GPU的所有线程，保证D矩阵中所有数据的一致性。

上述分解过程中同一个k值的情况下对L、U、D矩阵元素的计算过程可以使用GPU的多线程技术并行计算，每一个线程处理不同的(i,j)位置可以得到较好的并行效率。不同的k值则需要按照顺序计算，并使用线程同步操作同步共享内存中的数据。

本发明还提供了一种计算机可读存储介质，其上存储有计算机指令，所述计算机指令运行时执行上述的适用于GPU硬件的矩阵存储与计算方法的步骤，所述适用于GPU硬件的矩阵存储与计算方法参见前述部分的介绍，不再赘述。

本领域普通技术人员可以理解：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。