CN102567283B

CN102567283B - 一种利用gpu对小矩阵求逆的方法

Info

Publication number: CN102567283B
Application number: CN201110407357.8A
Authority: CN
Inventors: 隋丹; 李云洲; 周春晖; 赵熠飞; 赵明; 王京
Original assignee: Tsinghua University
Current assignee: Tsinghua University
Priority date: 2011-12-08
Filing date: 2011-12-08
Publication date: 2014-12-31
Anticipated expiration: 2031-12-08
Also published as: CN102567283A

Abstract

本发明公开了一种利用GPU对小矩阵求逆的方法，涉及无线通信领域。所述方法包括步骤：在GPU的共享存储器上创建维度为K×(N×N)的二维数组sm_a，创建两个维度为K×N的二维数组sm_is和sm_js；K和N均为大于0的自然数；将GPU的全局存储器中的K个N阶方阵并行存储到所述共享存储器的二维数组sm_a中；利用所述二维数组sm_is和sm_js，在所述共享存储器中完成对所述K个N阶方阵的求逆处理。所述方法既增加了并行的线程，又没有占用过多的共享存储器，且具有较好的可扩展性，显著提高了对小矩阵求逆运算的速度。

Description

一种利用GPU对小矩阵求逆的方法

技术领域

本发明涉及无线通信技术领域，特别涉及一种利用GPU对小矩阵求逆的方法。

背景技术

矩阵求逆是一种经常遇到的重要矩阵运算，在信号处理、神经网络、自动控制等领域都有广泛应用。特别是在4G无线通信标准中，多个关键功能模块，例如OFDM(Orthogonal Frequency DivisionMultiplexing，正交频分复用)系统信道估计、MIMO(Multiple-InputMultiple-Out-put，多输入输出天线系统)信号检测等，当采用迫零算法或最小均方误差算法时都可以归结为对信道矩阵的某种变换进行求逆操作，另外，对于码长较长的LDPC(Low Density Parity CheckCode，低密度奇偶校验码)码编码也需要进行大矩阵求逆。

矩阵求逆的处理速度直接影响了上述算法的执行速度，而矩阵求逆往往是非常费时的。现有的矩阵求逆大多是在CPU上通过软件实现的，能够满足较低数据传输速率的要求。也有部分在FPGA(Field-Programmable Gate Array，现场可编程门阵列)、DSP(Digital SignalProcessing，数字信号处理)硬件上实现矩阵求逆的，能够满足较高传输速率的要求，但灵活性、可配置性较差。近年来，随着GPU(Graphic Processing Unit，图形处理器)在非图形领域的科学计算中逐渐崭露头角，人们开始研究基于GPU的矩阵求逆算法。现有的基于GPU的矩阵求逆算法大多集中于高性能计算领域，针对维数较大(例如1024×1024)的矩阵，并且在一个应用中仅需进行一次大矩阵求逆。

而在无线通信系统中，需要对数量众多的小矩阵进行求逆处理。例如，LTE(Long Term Evolution，长期演进)标准中规定在带宽为5MHz时，可以采用2×2MIMO，或4×2MIMO，在20MHz带宽时，可以采用4×4MIMO，甚至8×8MIMO，此时信道矩阵的维度分别是2×2、4×2、4×4、8×8，经过变换后需要进行求逆处理的矩阵维度分别是2×2、4×4、8×8。而当带宽为5MHz、10MHz、15MHz、20MHz时，要求0.5ms的子帧周期内分别含有300、600、900、1200个OFDM符号，即要在0.5ms内分别完成300、600、900、1200个维度为2×2、4×4、8×8的矩阵的求逆处理。与对大矩阵的一次求逆相比，对大量小矩阵的求逆在算法流程、数据调度、计算线程和线程块的数据分发等方面都存在较大的不同。现有做法是，要么在一个计算线程中完成一个矩阵的求逆，要么在一个线程块中完成一个矩阵求逆。这两类做法相对比较直观，易于实现，但在GPU上的并行效率较低。这是因为，从GPU的硬件结构看，大量的CUDA(ComputeUnified Device Architecture，一种计算架构)核被分成若干个流多处理器(SMs)，例如最新的NVIDIA Tesla C2050由14个SM组成，每个SM包含32个CUDA核。每个SM作为一个单指令多线程(SIMT)处理器进行工作，而每个SM还含有一定大小的共享存储器，在共享存储器上的数据处理速度非常快，并且延时很小。而如果用一个线程计算一个矩阵的逆，那么每个线程上消耗的共享存储器较多，从而限制了SM上并发的线程个数，进而降低其并行效率。另一方面，如果用一个线程块计算一个矩阵的逆，即线程块中的每个线程处理矩阵的一个元素，由于我们要处理的矩阵尺寸往往较小(例如，2×2，4×4，8×8)，因此在一个线程块上的并行线程太小，也会影响其效率。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是：如何提供一种利用GPU对小矩阵求逆的方法，以提高对小矩阵求逆运算的速度。

(二)技术方案

为解决上述技术问题，本发明提供一种利用GPU对小矩阵求逆的方法，其包括步骤：

B：在GPU的共享存储器上创建维度为K×(N×N)的二维数组sm_a，创建两个维度为K×N的二维数组sm_is和sm_js；K和N均为大于0的自然数；

C：将GPU的全局存储器中的K个N阶方阵并行存储到所述共享存储器的二维数组sm_a中；

D：利用所述二维数组sm_is和sm_js，在所述共享存储器中完成对所述K个N阶方阵的求逆处理。

优选地，所述步骤D中，利用所述二维数组sm_is和sm_js，并采用全选主元高斯消去法，在所述共享存储器中并行完成对所述K个N阶方阵的求逆处理。

优选地，所述步骤D具体包括步骤：

D1：将K个N阶方阵A分别作为初始的当前方阵；

D2：判断K个当前方阵是否是1阶方阵，如果是，退出；否则，将K个当前方阵中的最大元素的行下标分别存储到所述二维数组sm_is中，列下标分别存储到所述二维数组sm_js中；

D3：对K个当前方阵，分别用所述行下标和列下标的组合对应的元素替换K个当前方阵中最上一行的对角线元素；

D4：对K个当前方阵中的非对角线元素根据如下公式按照从上至下并且从左至右的顺序进行更新：

A(k，j)＝A(k，j)/A(k，k)；

A(i，j)＝A(i，j)-A(i，k)×A(k，j)；

A(i，k)＝-A(i，k)/A(k，k)；

其中，0≤i，j≤N-1，i≠k，j≠k，i≠j；

D5：对K个当前方阵，分别删除最上一行和最左一列，得到新的K个当前方阵，执行步骤D2。

优选地，在所述步骤B之前还包括步骤A：选择由二维计算线程组成的线程块，所述线程块的第一维的数值对应待处理方阵的阶数，设定为N，第二维的数值对应待处理方阵的个数，设定为K。

优选地，在所述步骤D之后还包括步骤E：将所述K个N阶方阵的求逆结果从所述共享存储器转移到所述全局存储器。

优选地，所述N的取值为2、4或者8。

(三)有益效果

本发明所述利用GPU对小矩阵求逆的方法，让每个计算线程处理方阵的一行(或一列)的多个元素，一个线程块同时处理多个方阵，既增加了并行的线程，又没有占用过多的共享存储器，且具有较好的可扩展性，显著提高了对小矩阵求逆运算的速度。

附图说明

图1是本发明实施例所述的利用GPU对小矩阵求逆的方法流程图；

图2是本发明实施例所述的利用GPU对小矩阵求逆的方法加速效果图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

图1是本发明实施例所述的利用GPU对小矩阵求逆的方法流程图。如图1所示，所述方法包括：

步骤A：选择由二维计算线程组成的线程块，所述线程块的第一维的数值对应待处理方阵的阶数，设定为N，第二维的数值对应待处理方阵的个数，设定为K。N和K为大于0的自然数，其中，N的优选地取值为2、4或者8。

步骤B：在GPU的共享存储器上创建维度为K×(N×N)的二维数组sm_a，创建两个维度为K×N的二维数组sm_is和sm_js。

步骤C：将GPU的全局存储器中的K个N阶方阵A并行存储到所述共享存储器的二维数组sm_a中。所述二维数组sm_a可以按照行优先或者列优先的策略存储所述K个N阶方阵A中的元素。

步骤D：利用所述二维数组sm_is和sm_js，并采用全选主元高斯消去法，在所述共享存储器中并行完成对所述K个N阶方阵A的求逆处理。

所述步骤D具体包括：

步骤D1：将K个N阶方阵A分别作为初始的当前方阵。

步骤D2：判断K个当前方阵是否是1阶方阵，如果是，退出；否则，将K个当前方阵中的最大元素的行下标分别存储到所述二维数组sm_is中，列下标分别存储到所述二维数组sm_js中。

步骤D3：对K个当前方阵，分别用所述行下标和列下标的组合对应的元素替换K个当前方阵中最上一行的对角线元素。假设其中一个方阵A为4阶方阵，第一次循环时，其在二维数组sm_is中存储的行下标依次为1，其在二维数组sm_js中存储的列下标依次为2。则第一次执行该步骤D3中，用方阵A中的元素A(1，2)替换元素A(0，0)。

步骤D4：对K个当前方阵中的非对角线元素根据如下公式按照从上至下并且从左至右的顺序进行更新：

A(k，j)＝A(k，j)/A(k，k)；

A(i，j)＝A(i，j)-A(i，k)×A(k，j)；

A(i，k)＝-A(i，k)/A(k，k)；

其中，0≤i，j≤N-1，i≠k，j≠k，i≠j。

步骤D5：对K个当前方阵，分别删除最上一行和最左一列，得到新的K个当前方阵，执行步骤D2。

经过所述步骤D4和D5替换处理后得到的方阵是所述K个N阶方阵的求逆结果。

步骤E：将所述K个N阶方阵的求逆结果从所述共享存储器转移到所述全局存储器。

本实施例所述方法让每个计算线程处理方阵A的一行(或一列)的N个元素，一个线程块处理K个方阵，这样既增加了并行的线程，又没有占用过多的共享存储器。并且，所述方法能够灵活适用于不同阶次的方阵，具有较好的可扩展性。

为了测试加速结果，本实施例分别选取维度为2×2、4×4、8×8的方阵进行实验。实验中所采用的硬件配置如下：CPU为Intel Corei7-950(主频3.07GHz，内存6GB)；GPU为NVIDIATesla C2050(448个CUDA核处理器分为14个流多处理器，主频1.15GHz，显存3GB)；操作系统是Win764位专业版；编程环境为Visual Studio 2008；CUDA版本为4.0。为了便于描述加速结果，用T_CPU表示CPU进行矩阵求逆的运算时间，用T_GPU表示相应程序在GPU上的执行时间(包括GPU上核函数的运行时间及CPU和GPU之间数据拷贝时间的总和)，用T_CPU/T_GPU表示加速倍数。

表1是对比实验结果表。如表1所示，其给出了三种不同维度方阵进行10000次独立实验后统计得到的CPU与GPU运行时间比较结果，其中方阵数量均为60000。从表1中可以看出，对于三种维度的方阵，GPU的处理时间远远小于CPU的处理时间，并且方阵维度越小，加速比越高。

表1对比实验结果表

图2是本发明实施例所述的利用GPU对小矩阵求逆的方法加速效果图。该实验同样是针对维度为2×2、4×4、8×8的方阵，测试当方阵数量不同时，GPU对CPU的加速倍数。从图2中可以看出，对于同样个数的方阵而言，维度为2×2的方阵加速比最高，维度为8×8的方阵加速比最低。对于2×2的方阵，加速比随着处理方阵个数增加而快速提高，而对于维度为4×4、8×8的方阵，加速比受待处理方阵个数影响较小。

本发明实施例所述利用GPU对小矩阵求逆的方法，让每个计算线程处理方阵的一行(或一列)的多个元素，一个线程块同时处理多个方阵，既增加了并行的线程，又没有占用过多的共享存储器，且具有较好的可扩展性，显著提高了对小矩阵求逆运算的速度。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims

1.一种利用GPU对小矩阵求逆的方法，其特征在于，包括步骤：

A：选择由二维计算线程组成的线程块，所述线程块的第一维的数值对应待处理方阵的阶数，设定为N，第二维的数值对应待处理方阵的个数，设定为K；

C：将GPU的全局存储器中的K个N阶方阵并行存储到所述共享存储器的二维数组sm_a中；所述二维数组sm_a按照行优先或者列优先的策略存储所述K个N阶方阵中的元素；

D：利用所述二维数组sm_is和sm_js，在所述共享存储器中完成对所述K个N阶方阵的求逆处理；

所述步骤D具体包括步骤：

D1：将K个N阶方阵A分别作为初始的当前方阵；

A(k,j)＝A(k,j)/A(k,k)；

A(i,j)＝A(i,j)-A(i,k)×A(k,j)；

A(i,k)＝-A(i,k)/A(k,k)；

其中，0≤i,j≤N-1，i≠k，j≠k，i≠j；

D5：对K个当前方阵，分别删除最上一行和最左一列，得到新的K个当前方阵，执行步骤D2；

所述步骤D中，利用所述二维数组sm_is和sm_js，并采用全选主元高斯消去法，在所述共享存储器中并行完成对所述K个N阶方阵的求逆处理；

E：将所述K个N阶方阵的求逆结果从所述共享存储器转移到所述全局存储器；

所述N的取值为2、4或者8。