CN104572588A - 矩阵求逆处理方法和装置 - Google Patents

Abstract

本发明公开了一种矩阵求逆处理方法和装置。其中，矩阵求逆处理方法包括：获取经过中央处理器扩展后的扩展矩阵；根据扩展矩阵利用统一计算设备架构平台构建全局网格结构，其中，全局网格结构包括多个线程块及线程；利用全局网格结构对扩展矩阵的列向量进行并行处理，其中，对于扩展矩阵的每一列列向量上的数据采用串行计算的方式进行计算，得到计算结果，计算结果包括目标矩阵的逆矩阵和与逆矩阵具有相同大小的单位矩阵；以及将计算结果输出至中央处理器，其中，中央处理器从计算结果中提取逆矩阵。通过本发明，解决了现有技术中矩阵求逆过程中计算效率低的问题，达到了提高矩阵求逆过程中计算效率的效果。

Description

矩阵求逆处理方法和装置

技术领域

本发明涉及数据处理领域，具体而言，涉及一种矩阵求逆处理方法和装置。

背景技术

在现在数字化产业领域中，由于计算机技术的飞速发展，越来越多其他领域的技术开始被囊括进来，也越来越多的依靠数字化技术带来明显的产业推动能力，特别是制造业领域，由于设计、仿真分析在产品开发过程中所占有的比重增加，以及计算机便利的快速计算分析能力，现代制造业的产品开发流程几乎完全是基于数字化技术之上。也正是这种依赖关系，促使计算机软硬件不断更新换代以满足日益提升的性能需求。

对于矩阵的计算通常包括限单元法和有限差分法以及边界元法。相对于有限单元法和有限差分法，边界元法因其高精度及降维优势，更适合于快速预处理、自适应结构分析的过程。但采用边界元法得到的系数矩阵具有稠密且非对称的劣势，传统边界元法的求解效率以及求解问题规模都受到了限制。很多学者通过引入快速数值计算算法，如快速多极等，加速边界元求解计算过程，提高求解问题规模。然而，对于时间相关的工程问题，在时域或频域利用边界元法求解时，由于基本解的复杂性以及时域离散或转换需求，上述数值算法不能取得很好的计算效果。尤其是涉及O(N³)计算复杂度的矩阵求逆运算消耗了数值运算过程中的大量时间。

Hackbusch在1999年提出了基于遗传树结构的数值求解算法加速算法。利用矩阵分块，并将相容子块利用外乘法近似表示，从而压缩矩阵数据存储，减少参与矩阵运算的数据量，同时利用递归思想，借助遗传树结构，提出了近似求逆方法。但这种方法对于计算效率提升有限，同时由于是近似拟合，计算精度无法得到保证。此外，由于矩阵求逆运算中较强的数据关联性，难以利用中央处理器(CPU)多核并行计算提升降低计算时间耗费，这使得矩阵求逆的过程中计算效率较低。

针对现有技术中矩阵求逆过程中计算效率低的问题，目前尚未提出有效的解决方案。

发明内容

本发明的主要目的在于提供一种矩阵求逆处理方法和装置，以解决现有技术中矩阵求逆过程中计算效率低的问题。

为了实现上述目的，根据本发明实施例的一个方面，提供了一种矩阵求逆处理方法。根据本发明的矩阵求逆处理方法包括：获取经过中央处理器扩展后的扩展矩阵，所述扩展矩阵包括目标矩阵和与所述目标矩阵具有相同大小的单位矩阵；根据所述扩展矩阵利用统一计算设备架构平台构建全局网格结构，其中，所述全局网格结构包括多个线程；利用所述全局网格结构采用并行计算的方式依次对所述扩展矩阵的列向量进行计算，采用串行计算的方式对所述扩展矩阵的行向量的行间数据进行计算，得到计算结果，所述计算结果包括所述目标矩阵的逆矩阵和与所述逆矩阵具有相同大小的单位矩阵；以及将所述计算结果输出至所述中央处理器，其中，所述中央处理器从所述计算结果中提取所述逆矩阵。

进一步地，根据所述扩展矩阵利用统一计算设备架构平台构建全局网格结构包括：根据所述扩展矩阵的大小确定分割基数；按照所述分割基数对所述扩展矩阵的行向量和列向量进行划分，得到多个数据段；以及根据所述分割基数和所述多个数据段的数量构建线程块结构，形成所述全局网格结构，其中，所述全局网格结构包括与所述多个数据段一一对应的线程块，每个线程块具有与所述分割基数相同数量的线程。

进一步地，利用所述全局网格结构采用并行计算的方式依次对所述扩展矩阵的列向量进行计算，采用串行计算的方式对所述扩展矩阵的行向量的行间数据进行计算，得到计算结果包括：计算所述扩展矩阵中当前行向量的系数向量，所述系数向量包括所述扩展矩阵中所述当前行向量之外的其他行向量对于所述当前行向量的系数；利用所述全局网格结构中线程计算该线程在所述扩展矩阵的映射位置上的变换结果，用所述变换结果替换在所述扩展矩阵的映射位置上的数据，得到替换后的扩展矩阵；判断所述当前行向量是否是所述扩展矩阵的最后一行行向量；如果判断出所述当前行向量不是所述扩展矩阵的最后一行行向量，则将所述当前行向量的下一行行向量作为所述当前行向量，返回执行计算当前行向量的系数向量的步骤；以及如果判断出所述当前行向量是所述扩展矩阵的最后一行行向量，则对所述替换后的扩展矩阵进行单位化处理，得到所述计算结果。

进一步地，计算所述扩展矩阵中当前行向量的系数向量包括：获取所述当前行向量的对角位置数据；获取所述对角位置数据所在的列向量上的数据；依次将对角位置数据与所述对角位置数据所在的列向量上的数据相除，得到所述系数向量。

进一步地，在利用所述全局网格结构中线程计算该线程在所述扩展矩阵的映射位置上的变换结果之前，所述矩阵求逆处理方法还包括：建立共享显存空间，所述共享显存空间用于存储所述当前行向量的数据；其中，所述全局网格结构中线程从所述共享显存空间调用所述当前行向量的数据用于计算该线程在所述扩展矩阵的映射位置上的变换结果。

为了实现上述目的，根据本发明实施例的另一方面，提供了一种矩阵求逆处理装置。根据本发明的矩阵求逆处理装置包括：获取单元，用于获取经过中央处理器扩展后的扩展矩阵，所述扩展矩阵包括目标矩阵和与所述目标矩阵具有相同大小的单位矩阵；第一建立单元，用于根据所述扩展矩阵利用统一计算设备架构平台构建全局网格结构，其中，所述全局网格结构包括多个线程；计算单元，用于利用所述全局网格结构采用并行计算的方式依次对所述扩展矩阵的列向量进行计算，采用串行计算的方式对所述扩展矩阵的行向量的行间数据进行计算，得到计算结果，所述计算结果包括所述目标矩阵的逆矩阵和与所述逆矩阵具有相同大小的单位矩阵；以及输出单元，用于将所述计算结果输出至所述中央处理器，其中，所述中央处理器从所述计算结果中提取所述逆矩阵。

进一步地，所述第一建立单元包括：确定模块，用于根据所述扩展矩阵的大小确定分割基数；划分模块，用于按照所述分割基数对所述扩展矩阵的行向量和列向量进行划分，得到多个数据段；以及构建模块，用于根据所述分割基数和所述多个数据段的数量构建线程块结构，形成所述全局网格结构，其中，所述全局网格结构包括与所述多个数据段一一对应的线程块，每个线程块具有与所述分割基数相同数量的线程。

进一步地，所述计算单元包括：第一计算模块，用于计算所述扩展矩阵中当前行向量的系数向量，所述系数向量包括所述扩展矩阵中所述当前行向量之外的其他行向量对于所述当前行向量的系数；第二计算模块，用于利用所述全局网格结构中线程计算该线程在所述扩展矩阵的映射位置上的变换结果，用所述变换结果替换在所述扩展矩阵的映射位置上的数据，得到替换后的扩展矩阵；判断模块，用于判断所述当前行向量是否是所述扩展矩阵的最后一行行向量；所述第一计算模块还用于如果判断出所述当前行向量不是所述扩展矩阵的最后一行行向量，则将所述当前行向量的下一行行向量作为所述当前行向量，计算当前行向量的系数向量；以及处理模块，用于如果判断出所述当前行向量是所述扩展矩阵的最后一行行向量，则对所述替换后的扩展矩阵进行单位化处理，得到所述计算结果。

进一步地，所述第一计算模块包括：第一获取子模块，用于获取所述当前行向量的对角位置数据；第二获取子模块，用于获取所述对角位置数据所在的列向量上的数据；除法子模块，用于依次将对角位置数据与所述对角位置数据所在的列向量上的数据相除，得到所述系数向量。

进一步地，所述矩阵求逆处理装置还包括：第二建立单元，用于在利用所述全局网格结构中线程计算该线程在所述扩展矩阵的映射位置上的变换结果之前，建立共享显存空间，所述共享显存空间用于存储所述当前行向量的数据；所述第二计算模块还用于所述全局网格结构中线程从所述共享显存空间调用所述当前行向量的数据用于计算该线程在所述扩展矩阵的映射位置上的变换结果。

根据本发明实施例，获取经过中央处理器扩展后的扩展矩阵，根据扩展矩阵利用统一计算设备架构平台构建全局网格结构，利用全局网格结构采用并行计算的方式依次对扩展矩阵的列向量进行计算，采用串行计算的方式对扩展矩阵的行向量的行间数据进行计算，得到计算结果，将计算结果输出至中央处理器，其中，中央处理器从计算结果中提取逆矩阵，解决了现有技术中矩阵求逆过程中计算效率低的问题，达到了提高矩阵求逆过程中计算效率的效果。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是根据本发明实施例的矩阵求逆处理方法的流程图；

图2是根据本发明实施例的全局网格结构的示意图；

图3是根据本发明实施例优选的矩阵求逆处理方法的流程图；以及

图4是根据本发明实施例的矩阵求逆处理装置的示意图。

具体实施方式

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本发明。

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

需要说明的是，本发明的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本发明的实施例。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、***、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

本发明实施例提供了一种矩阵求逆处理方法。该方法可以由图形处理器(GPU)来执行。

图1是根据本发明实施例的矩阵求逆处理方法的流程图。如图1所示，该矩阵求逆处理方法包括步骤如下：

步骤S102，获取经过中央处理器扩展后的扩展矩阵，扩展矩阵包括目标矩阵和与目标矩阵具有相同大小的单位矩阵。

目标矩阵即为待求逆矩阵。在计算目标矩阵之前，先通过主机端中央处理器(CPU)将目标矩阵扩展为扩展矩阵，例如m*n大小的目标矩阵[M]，将其扩展为扩展矩阵[M|E]，然后将该扩展矩阵输出至设备端图形处理器(GPU)，通过GPU来实现目标矩阵求逆。

步骤S104，根据扩展矩阵利用统一计算设备架构平台构建全局网格结构。其中，全局网格结构包括多个线程。

统一计算设备架构(CUDA)平台是C语言的一种扩展，它允许使用标准C来进行GPU代码编程。这个代码既适用于中央处理器(CPU)，也适用于图形处理器(GPU)。主机端负责派生出运行在GPU设备端上的多线程任务(kernel函数)；设备端设有内部调度器来将kernel程序分配到相应的GPU硬件上。利用CUDA平台构建由多个线程全组成的局网格结构，利用该全局网格结构对扩展矩阵进行处理。具体地，全局网格结构可以包括至少三个子全局网格结构，分别用于系数向量计算、矩阵行运算和矩阵单位化处理三个过程，其中，系数向量计算对应的子全局网格结构在横、纵方向划分线程块，其中每一个线程用于处理扩展矩阵中目标矩阵上的每一个数据元素，矩阵行运算对应的子全局网格结构仅在横向划分线程块，其中每一个线程用于处理扩展矩阵的一个列向量，矩阵单位化处理对应的子全局网格结构仅在纵向划分线程块，其中每一个线程用于处理扩展矩阵的一个行向量。

步骤S106，利用所述全局网格结构对所述扩展矩阵的列向量进行并行处理，其中，对于所述扩展矩阵的每一列列向量上的数据采用串行计算的方式进行计算，得到计算结果。计算结果包括目标矩阵的逆矩阵和与逆矩阵具有相同大小的单位矩阵。

由于在进行目标矩阵求逆的过程中，扩展矩阵的行向量与行向量之间的数据即行间数据的关联性较高，因此，对扩展矩阵的不同行向量之间的数据采用串行计算的方式进行计算，即在每个线程内按照行向量的顺序串行计算行间数据。而对于扩展矩阵的每个列向量，可以采用并行计算，即每个线程之间并行计算，从而提高行间数据的计算效率。

步骤S108，将计算结果输出至中央处理器，其中，中央处理器从计算结果中提取逆矩阵。

设备端GPU计算得到计算结果之后，将该计算结果输出给中央处理器CPU，通过CPU提取得到目标矩阵的逆矩阵。

以m*n大小的目标矩阵[M]为例，在计算目标矩阵的逆矩阵之前，需要根据目标矩阵[M]大小分配数据存储空间。对于主机端，需要分配m*n大小空间存储目标矩阵；在设备端，则需要分配2*m*n*sizeof(float)大小空间(g_iMatrix)存储中间计算数据，同时开辟sizeof(float)*m*1大小空间(g_tVector)存储计算向量。然后将主机端的目标矩阵数据传入显存空间，利用GPU对目标矩阵进行求逆。

CUDA显存中求得的矩阵实际为[E|InM]，其中InM为原矩阵的逆阵，[E]是m*n大小的单位矩阵。将数据由显存空间传入内存空间，并提取InM替换原矩阵M，同时将显存及内存空间释放。

根据本发明实施例，获取经过中央处理器扩展后的扩展矩阵，根据扩展矩阵利用统一计算设备架构平台构建全局网格结构，利用所述全局网格结构对所述扩展矩阵的列向量进行并行处理，其中，对于所述扩展矩阵的每一列列向量上的数据采用串行计算的方式进行计算，得到计算结果，将计算结果输出至中央处理器，其中，中央处理器从计算结果中提取逆矩阵。利用CPU主机端强大的逻辑处理能力和GPU设备端强大的运算能力，通过将矩阵求逆过程的行、列数据处理过程分拆，最大化利用CUDA的并行计算能力降低计算时间，提升矩阵求逆计算效率。

优选地，根据扩展矩阵利用统一计算设备架构平台构建全局网格结构包括：根据扩展矩阵的大小确定分割基数；按照分割基数对扩展矩阵的行向量和列向量进行划分，得到多个数据段；以及根据分割基数和多个数据段的数量构建线程块结构，形成全局网格结构，其中，全局网格结构包括与多个数据段一一对应的线程块，每个线程块具有与分割基数相同数量的线程。

具体地，全局网格结构可以包括至少三个子全局网格结构，分别用于系数向量计算、矩阵行运算和矩阵单位化处理三个过程。

以m*n大小的目标矩阵[M]为例，扩展为2*m*n大小的矩阵[M|E]。

将行m及列n按照分割基数s(可以取s＝8、16、32或64，按照具体计算数据规模大小选择)对应分为Segment_m＝m/s+1、Segment_n＝n/s+1大小的数据段。利用行和列分割得到的Segments构建全局Grid(网格)；根据分割基数s，构建Block结构，形成全局网格结构。如图2所示，T表示线程(Thread)，B表示线程块(Block)，x，y表示线程块在两个方向上的维度，X，Y表示网格两个方向上的维度。

具体地，计算过程中，在X，Y，Z三个维度的构建原则为：系数计算时的Grid结构为(Segment_m,Segment_n,1)，Block结构为(s,s,1)；矩阵行运算的Grid结构为(1,2*Segment_n,1)，Block结构为(1,s,1)；矩阵单位化的Grid结构为Grid(Segment_m,1,1)，Block结构为(s,1,1)。

优选地，利用所述全局网格结构对所述扩展矩阵的列向量进行并行处理，其中，对于所述扩展矩阵的每一列列向量上的数据采用串行计算的方式进行计算，得到计算结果包括：计算扩展矩阵中当前行向量的系数向量，系数向量包括扩展矩阵中当前行向量之外的其他行向量对于当前行向量的系数；利用全局网格结构中线程计算该线程在扩展矩阵的映射位置上的变换结果，用变换结果替换在扩展矩阵的映射位置上的数据，得到替换后的扩展矩阵；判断当前行向量是否是扩展矩阵的最后一行行向量；如果判断出当前行向量不是扩展矩阵的最后一行行向量，则将当前行向量的下一行行向量作为当前行向量，返回执行计算当前行向量的系数向量的步骤；以及如果判断出当前行向量是扩展矩阵的最后一行行向量，则对替换后的扩展矩阵进行单位化处理，得到计算结果。

具体地，可以通过串行过程遍历扩展矩阵行向量k(k行为扩展矩阵中遍历到的扩展矩阵的当前行向量)，借助CUDA线程对扩展矩阵的列进行并行运算，计算其他行对于k行的系数，得到系数向量，存入向量空间g_tVector。

针对当前遍历的k行，利用系数向量，在每个线程中计算矩阵映射位置中的变换结果，即是矩阵行运算，并替换原数据。

优选地，在利用全局网格结构中线程计算该线程在扩展矩阵的映射位置上的变换结果之前，矩阵求逆处理方法还包括：建立共享显存空间，共享显存空间用于存储当前行向量的数据，其中，全局网格结构中线程从共享显存空间调用当前行向量的数据用于计算该线程在扩展矩阵的映射位置上的变换结果。

通过建立共享显存空间来存储当前行向量的数据，以便于在计算过程中，不同线程可以同时调用该行向量的数据，避免重复调用占用计算资源。

在每个线程中计算原矩阵数据和当前矩阵对角位置数据(base)进行除法运算，并替换原数据，即进行矩阵映射位置数据的单位化操作。最后，判断是否遍历到扩展矩阵的最后一行，否则继续对下一行行向量执行上述步骤。

根据本发明实施例，通过依次遍历矩阵行向量，各行向量的数据采用并行计算，从而提高矩阵求逆的计算效率。

进一步地，计算扩展矩阵中当前行向量的系数向量包括：获取当前行向量的对角位置数据；获取对角位置数据所在的列向量上的数据；依次将对角位置数据与对角位置数据所在的列向量上的数据相除，得到系数向量。

具体地，在系数计算时的Grid结构为(Segment_m,Segment_n,1)，Block结构为(s,s,1)。在Grid和Block结构下，每个线程处理全局矩阵(g_iMatrix)相应位置数据。在线程核函数中，线程映射在全局中的X、Y或Z方向位置(tid_in_grid_x/y/z)由Block位置与宽度和Block中线程位置确定，如:

tid_in_grid_y＝blockDim.y*blockIdx.y+threadIdx.y。

查找当前线程在Y方向位置所映射的矩阵所在行，并将此行的对角位置数据与tid_in_grid_y行所在的矩阵对角位置数据(base)相除，计算结果放入系数向量(g_tVector)对应tid_in_grid_y的位置中。执行伪代码如下：

在矩阵行运算时，Grid结构为(1,Segment_m,1)，Block结构为(1,s,1)。

利用线程全局位置确定线程所要处理的矩阵数据位置(tid_in_grid)，同时确定需要进行减法计算的k行对应列位置(T_in_obj)。

tid_in_grid＝tid_in_grid_y*n+tid_in_grid_x；

T_in_obj＝k*n+tid_in_grid_x；

由于k行数据是每个线程计算中均需要调用的数据，为避免重复调用占用计算资源，在线程块中建立共享显存空间(sdata)存储这些数据，用于不同线程同时访问。

矩阵的k行、T_in_obj列的数据和g_iVector中对应位置的系数相乘，并与tid_in_grid位置对应的矩阵数据进行减法运算，并进行数据替换。

在矩阵单位化计算过程中，Grid结构为Grid(Segment_m,1,1)，Block结构为(s,1,1)。在任一线程中，索引第k行并对应tid_in_grid_x列的矩阵位置tid_in_matrix，将此位置数据除以当前矩阵base对应在k行位置的数据，将结果存入矩阵tid_in_matrix位置。

以大小为100*100目标矩阵[M]为例，结合图3来对本发明实施例进行描述。

步骤S302，根据目标矩阵[M]大小分配数据存储空间。对于主机端，需要分配100*100大小空间存储目标矩阵；在设备端，则需要分配156.25KB大小空间(g_iMatrix)存储中间计算数据，同时开辟1KB大小空间(g_tVector)存储计算向量。将主机端待求逆矩阵数据传入显存空间。

步骤S304，将目标矩阵[M]扩展为矩阵[M|E]。其中矩阵[M|E]大小为2*100*100。

步骤S306，将行m及列n分为相同大小的若干数据段。具体地，将行m及列n按照分割基数s(s＝8、16、32或64，按照具体计算数据规模小选择)对应分为Segment_m＝7、Segment_n＝7大小的数据段。

步骤S308，利用行和列分割得到的数据段构建全局网格结构。具体地，根据分割基数s，构建Block结构。具体在X,Y,Z三个维度的构建原则为：系数计算时的Grid结构为(7,7,1)，Block结构为(16,16,1)；矩阵行运算的Grid结构为(1,14,1)，Block结构为(1,16,1)；矩阵单位化的Grid结构为Grid(7,1,1)，Block结构为(16,1,1)。

在构建全局网格结构之后，在每个线程中并行执行以下步骤S310至步骤S316，即对扩展矩阵的列向量并行计算。

步骤S310，遍历矩阵的当前行向量k。具体地，从第1行至第100行，循环遍历矩阵行向量，k＝1,2,…100。

步骤S312，系数计算。在Grid和Block结构下，每个线程处理全局矩阵(g_iMatrix)相应位置数据。在线程核函数中，线程映射在全局中的X、Y或Z方向位置(tid_in_grid_x/y/z)由Block位置与宽度和Block中线程位置确定，如:

tid_in_grid_y＝blockDim.y*blockIdx.y+threadIdx.y。

查找当前线程在Y方向位置所映射的矩阵所在行，并将此行的对角位置数据与tid_in_grid_y行所在的矩阵对角位置数据(base)相除，计算结果放入系数向量(g_tVector)对应tid_in_grid_y的位置中。执行伪代码如下：

步骤S314，矩阵行运算，CUDA中每一线程数据和当前遍历的矩阵k行数据进行运算。首先利用线程全局位置确定线程所要处理的矩阵数据位置(tid_in_grid)，同时确定需要进行减法计算的k行对应列位置(T_in_obj)。

tid_in_grid＝tid_in_grid_y*n+tid_in_grid_x；

T_in_obj＝k*n+tid_in_grid_x；

由于k行数据是每个线程计算中均需要调用的数据，为避免重复调用占用计算资源，在线程块中建立共享显存空间(sdata)存储这些数据，用于不同线程同时访问。

矩阵的k行、T_in_obj列的数据和g_iVector中对应位置的系数相乘，并与tid_in_grid位置对应的矩阵数据进行减法运算，并进行数据替换。

步骤S316，矩阵单位化。在任一线程中，索引第k行并对应tid_in_grid_x列的矩阵位置tid_in_matrix，将此位置数据除以当前矩阵base对应在k行位置的数据，将计算结果存入矩阵tid_in_matrix位置。

步骤S318，判断k是否等于m。如果是，则k加1，并返回步骤S310，反之，则执行步骤S320。

步骤S320，CUDA显存中求得的矩阵实际为[E|InM]，其中InM为原矩阵的逆阵，[E]是100*100大小的单位矩阵。将数据由显存空间传入内存空间，并提取InM替换原矩阵M，同时将显存及内存空间释放。

本发明实施例还提供了一种矩阵求逆处理装置。需要说明的是，本发明实施例的矩阵求逆处理装置可以用于执行本发明实施例所提供的矩阵求逆处理方法，本发明实施例的矩阵求逆处理方法也可以通过本发明实施例所提供的矩阵求逆处理装置来执行。

图4是根据本发明实施例的矩阵求逆处理装置的示意图。如图4所示，该矩阵求逆处理装置包括：获取单元10、第一建立单元20、计算单元30和输出单元40。

获取单元10用于获取经过中央处理器扩展后的扩展矩阵，扩展矩阵包括目标矩阵和与目标矩阵具有相同大小的单位矩阵。

目标矩阵即为待求逆矩阵。在计算目标矩阵之前，先通过主机端中央处理器(CPU)将目标矩阵扩展为扩展矩阵，例如m*n大小的目标矩阵[M]，将其扩展为扩展矩阵[M|E]，然后将该扩展矩阵输出至设备端图形处理器(GPU)，通过GPU来实现目标矩阵求逆。

第一建立单元20用于根据扩展矩阵利用统一计算设备架构平台构建全局网格结构。其中，全局网格结构包括多个线程。

统一计算设备架构(CUDA)平台，CUDA是C语言的一种扩展，它允许使用标准C来进行GPU代码编程。这个代码既适用于中央处理器(CPU)，也适用于图形处理器(GPU)。主机端负责派生出运行在GPU设备端上的多线程任务(kernel函数)；设备端设有内部调度器来将kernel程序分配到相应的GPU硬件上。利用CUDA平台构建由多个线程全组成的局网格结构，利用该全局网格结构对扩展矩阵进行处理。具体地，全局网格结构可以包括至少三个子全局网格结构，分别用于系数向量计算、矩阵行运算和矩阵单位化处理三个过程，其中，系数向量计算对应的子全局网格结构在横、纵方向划分线程块，其中每一个线程用于处理扩展矩阵中目标矩阵上的每一个数据元素，矩阵行运算对应的子全局网格结构仅在横向划分线程块，其中每一个线程用于处理扩展矩阵的一个列向量，矩阵单位化处理对应的子全局网格结构仅在纵向划分线程块，其中每一个线程用于处理扩展矩阵的一个行向量。

计算单元30用于利用所述全局网格结构对所述扩展矩阵的列向量进行并行处理，其中，对于所述扩展矩阵的每一列列向量上的数据采用串行计算的方式进行计算，得到计算结果。计算结果包括目标矩阵的逆矩阵和与逆矩阵具有相同大小的单位矩阵。

由于在进行目标矩阵求逆的过程中，扩展矩阵的行向量与行向量之间的数据即行间数据的关联性较高，因此，对扩展矩阵的不同行向量之间的数据采用串行计算的方式进行计算，即在每个线程内按照行向量的顺序串行计算行间数据。而对于扩展矩阵的每个列向量，可以采用并行计算，即每个线程之间并行计算，从而提高行间数据的计算效率。

输出单元40用于将计算结果输出至中央处理器，其中，中央处理器从计算结果中提取逆矩阵。

设备端GPU计算得到计算结果之后，将该计算结果输出给中央处理器CPU，通过CPU提取得到目标矩阵的逆矩阵。

以m*n大小的目标矩阵[M]为例，在计算目标矩阵的逆矩阵之前，需要根据目标矩阵[M]大小分配数据存储空间。对于主机端，需要分配m*n大小空间存储目标矩阵；在设备端，则需要分配2*m*n*sizeof(float)大小空间(g_iMatrix)存储中间计算数据，同时开辟sizeof(float)*m*1大小空间(g_tVector)存储计算向量。然后将主机端的目标矩阵数据传入显存空间，利用GPU对目标矩阵进行求逆。

CUDA显存中求得的矩阵实际为[E|InM]，其中InM为原矩阵的逆阵，[E]是m*n大小的单位矩阵。将数据由显存空间传入内存空间，并提取InM替换原矩阵M，同时将显存及内存空间释放。

根据本发明实施例，获取经过中央处理器扩展后的扩展矩阵，根据扩展矩阵利用统一计算设备架构平台构建全局网格结构，利用所述全局网格结构对所述扩展矩阵的列向量进行并行处理，其中，对于所述扩展矩阵的每一列列向量上的数据采用串行计算的方式进行计算，得到计算结果，将计算结果输出至中央处理器，其中，中央处理器从计算结果中提取逆矩阵。利用CPU主机端强大的逻辑处理能力和GPU设备端强大的运算能力，通过将矩阵求逆过程的行、列数据处理过程分拆，最大化利用CUDA的并行计算能力降低计算时间，提升矩阵求逆计算效率。

优选地，第一建立单元包括：确定模块，用于根据扩展矩阵的大小确定分割基数；划分模块，用于按照分割基数对扩展矩阵的行向量和列向量进行划分，得到多个数据段；以及构建模块，用于根据分割基数和多个数据段的数量构建线程块结构，形成全局网格结构，其中，全局网格结构包括与多个数据段一一对应的线程块，每个线程块具有与分割基数相同数量的线程。

具体地，全局网格结构可以包括至少三个子全局网格结构，分别用于系数向量计算、矩阵行运算和矩阵单位化处理三个过程。

以m*n大小的目标矩阵[M]为例，扩展为2*m*n大小的矩阵[M|E]。

将行m及列n按照分割基数s(可以取s＝8、16、32或64，按照具体计算数据规模大小选择)对应分为Segment_m＝m/s+1、Segment_n＝n/s+1大小的数据段。利用行和列分割得到的Segments构建全局Grid(网格)；根据分割基数s，构建Block结构，形成全局网格结构。如图2所示，T表示线程(Thread)，B表示线程块(Block)，x，y表示线程块在两个方向上的维度，X，Y表示网格两个方向上的维度。

具体地，计算过程中，在X，Y，Z三个维度的构建原则为：系数计算时的Grid结构为(Segment_m,Segment_n,1)，Block结构为(s,s,1)；矩阵行运算的Grid结构为(1,2*Segment_n,1)，Block结构为(1,s,1)；矩阵单位化的Grid结构为Grid(Segment_m,1,1)，Block结构为(s,1,1)。

优选地，计算单元包括：第一计算模块，用于计算扩展矩阵中当前行向量的系数向量，系数向量包括扩展矩阵中当前行向量之外的其他行向量对于当前行向量的系数；第二计算模块，用于利用全局网格结构中线程计算该线程在扩展矩阵的映射位置上的变换结果，用变换结果替换在扩展矩阵的映射位置上的数据，得到替换后的扩展矩阵；判断模块，用于判断当前行向量是否是扩展矩阵的最后一行行向量；第一计算模块还用于如果判断出当前行向量不是扩展矩阵的最后一行行向量，则将当前行向量的下一行行向量作为当前行向量，计算当前行向量的系数向量；以及处理模块，用于如果判断出当前行向量是扩展矩阵的最后一行行向量，则对替换后的扩展矩阵进行单位化处理，得到计算结果。

具体地，可以通过串行过程遍历扩展矩阵行向量k(k行为扩展矩阵中遍历到的扩展矩阵的当前行向量)，借助CUDA线程对扩展矩阵的列进行并行运算，计算其他行对于k行的系数，得到系数向量，存入向量空间g_tVector。

针对当前遍历的k行，利用系数向量，在每个线程中计算矩阵映射位置中的变换结果，即是矩阵行运算，并替换原数据。

优选地，矩阵求逆处理装置还包括：第二建立单元，用于在利用全局网格结构中线程计算该线程在扩展矩阵的映射位置上的变换结果之前，建立共享显存空间，共享显存空间用于存储当前行向量的数据；第二计算模块还用于全局网格结构中线程从共享显存空间调用当前行向量的数据用于计算该线程在扩展矩阵的映射位置上的变换结果。

通过建立共享显存空间来存储当前行向量的数据，以便于在计算过程中，不同线程可以同时调用该行向量的数据，避免重复调用占用计算资源。

在每个线程中计算原矩阵数据和当前矩阵对角位置数据(base)进行除法运算，并替换原数据，即进行矩阵映射位置数据的单位化操作。最后，判断是否遍历到扩展矩阵的最后一行，否则继续对下一行行向量执行上述步骤。

根据本发明实施例，通过依次遍历矩阵行向量，各行向量的数据采用并行计算，从而提高矩阵求逆的计算效率。

进一步地，第一计算模块包括：第一获取子模块，用于获取当前行向量的对角位置数据；第二获取子模块，用于获取对角位置数据所在的列向量上的数据；除法子模块，用于依次将对角位置数据与对角位置数据所在的列向量上的数据相除，得到系数向量。

具体地，在系数计算时的Grid结构为(Segment_m,Segment_n,1)，Block结构为(s,s,1)。在Grid和Block结构下，每个线程处理全局矩阵(g_iMatrix)相应位置数据。在线程核函数中，线程映射在全局中的X、Y或Z方向位置(tid_in_grid_x/y/z)由Block位置与宽度和Block中线程位置确定，如:

tid_in_grid_y＝blockDim.y*blockIdx.y+threadIdx.y。

查找当前线程在Y方向位置所映射的矩阵所在行，并将此行的对角位置数据与tid_in_grid_y行所在的矩阵对角位置数据(base)相除，计算结果放入系数向量(g_tVector)对应tid_in_grid_y的位置中。执行伪代码如下：

在矩阵行运算时，Grid结构为(1,Segment_m,1)，Block结构为(1,s,1)。

利用线程全局位置确定线程所要处理的矩阵数据位置(tid_in_grid)，同时确定需要进行减法计算的k行对应列位置(T_in_obj)。

tid_in_grid＝tid_in_grid_y*n+tid_in_grid_x；

T_in_obj＝k*n+tid_in_grid_x；

由于k行数据是每个线程计算中均需要调用的数据，为避免重复调用占用计算资源，在线程块中建立共享显存空间(sdata)存储这些数据，用于不同线程同时访问。

矩阵的k行、T_in_obj列的数据和g_iVector中对应位置的系数相乘，并与tid_in_grid位置对应的矩阵数据进行减法运算，并进行数据替换。

在矩阵单位化计算过程中，Grid结构为Grid(Segment_m,1,1)，Block结构为(s,1,1)。在任一线程中，索引第k行并对应tid_in_grid_x列的矩阵位置tid_in_matrix，将此位置数据除以当前矩阵base对应在k行位置的数据，将结果存入矩阵tid_in_matrix位置。

本发明的效果可以通过以下仿真和实测数据实验进一步说明。

仿真条件

算法运行平台:

CPU:Intel(R)Xeon(R)CPU E5-1620v2(3.70GHz)；

GPU:NVIDIA Quadro K4000；

内存:16GB；

编译器:Visual Studio 2010；

仿真内容

针对边界单元法数值计算中的稠密非对称矩阵分别进行串行及CUDA并行的求逆运算，统计计算花费时间(单位为秒)，并验证计算结果精度。

实测数据实验

求解了两个稠密非堆成矩阵，规模分别是3084*3084以及7605*7605。对于第一种求解规模，并行CUDA求逆花费时间为：5.383s；串行CPU求逆花费时间为：32.325s；加速比为6.01。

对于第一种求解规模，并行CUDA求逆花费时间为：14.447；串行CPU求逆花费时间为：198.144；加速比为13.72。

可以看到，相对于串行CPU计算，CUDA并行求逆算法的计算效率提升明显，并随着计算数据规模提升而提高。

需要说明的是，对于前述的各方法实施例，为了简单描述，故将其都表述为一系列的动作组合，但是本领域技术人员应该知悉，本发明并不受所描述的动作顺序的限制，因为依据本发明，某些步骤可以采用其他顺序或者同时进行。其次，本领域技术人员也应该知悉，说明书中所描述的实施例均属于优选实施例，所涉及的动作和模块并不一定是本发明所必须的。

在上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详述的部分，可以参见其他实施例的相关描述。

在本申请所提供的几个实施例中，应该理解到，所揭露的装置，可通过其它的方式实现。例如，以上所描述的装置实施例仅仅是示意性的，例如所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个***，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，装置或单元的间接耦合或通信连接，可以是电性或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。

所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可为个人计算机、移动终端、服务器或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、移动硬盘、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种矩阵求逆处理方法，其特征在于，包括：

获取经过中央处理器扩展后的扩展矩阵，所述扩展矩阵包括目标矩阵和与所述目标矩阵具有相同大小的单位矩阵；

根据所述扩展矩阵利用统一计算设备架构平台构建全局网格结构，其中，所述全局网格结构包括多个线程；

利用所述全局网格结构对所述扩展矩阵的列向量进行并行处理，其中，对于所述扩展矩阵的每一列列向量上的数据采用串行计算的方式进行计算，得到计算结果，所述计算结果包括所述目标矩阵的逆矩阵和与所述逆矩阵具有相同大小的单位矩阵；以及

将所述计算结果输出至所述中央处理器，其中，所述中央处理器从所述计算结果中提取所述逆矩阵。

2.根据权利要求1所述的矩阵求逆处理方法，其特征在于，根据所述扩展矩阵利用统一计算设备架构平台构建全局网格结构包括：

根据所述扩展矩阵的大小确定分割基数；

按照所述分割基数对所述扩展矩阵的行向量和列向量进行划分，得到多个数据段；以及

根据所述分割基数和所述多个数据段的数量构建线程块结构，形成所述全局网格结构，其中，所述全局网格结构包括与所述多个数据段一一对应的线程块，每个线程块具有与所述分割基数相同数量的线程。

3.根据权利要求1所述的矩阵求逆处理方法，其特征在于，利用所述全局网格结构对所述扩展矩阵的列向量进行并行处理，其中，对于所述扩展矩阵的每一列列向量上的数据采用串行计算的方式进行计算，得到计算结果包括：

计算所述扩展矩阵中当前行向量的系数向量，所述系数向量包括所述扩展矩阵中所述当前行向量之外的其他行向量对于所述当前行向量的系数；

利用所述全局网格结构中线程计算该线程在所述扩展矩阵的映射位置上的变换结果，用所述变换结果替换在所述扩展矩阵的映射位置上的数据，得到替换后的扩展矩阵；

判断所述当前行向量是否是所述扩展矩阵的最后一行行向量；

如果判断出所述当前行向量不是所述扩展矩阵的最后一行行向量，则将所述当前行向量的下一行行向量作为所述当前行向量，返回执行计算当前行向量的系数向量的步骤；以及

如果判断出所述当前行向量是所述扩展矩阵的最后一行行向量，则对所述替换后的扩展矩阵进行单位化处理，得到所述计算结果。

4.根据权利要求3所述的矩阵求逆处理方法，其特征在于，计算所述扩展矩阵中当前行向量的系数向量包括：

获取所述当前行向量的对角位置数据；

获取所述对角位置数据所在的列向量上的数据；

依次将对角位置数据与所述对角位置数据所在的列向量上的数据相除，得到所述系数向量。

5.根据权利要求3所述的矩阵求逆处理方法，其特征在于，在利用所述全局网格结构中线程计算该线程在所述扩展矩阵的映射位置上的变换结果之前，所述矩阵求逆处理方法还包括：

建立共享显存空间，所述共享显存空间用于存储所述当前行向量的数据；

其中，所述全局网格结构中线程从所述共享显存空间调用所述当前行向量的数据用于计算该线程在所述扩展矩阵的映射位置上的变换结果。

6.一种矩阵求逆处理装置，其特征在于，包括：

获取单元，用于获取经过中央处理器扩展后的扩展矩阵，所述扩展矩阵包括目标矩阵和与所述目标矩阵具有相同大小的单位矩阵；

第一建立单元，用于根据所述扩展矩阵利用统一计算设备架构平台构建全局网格结构，其中，所述全局网格结构包括多个线程；

计算单元，用于利用所述全局网格结构对所述扩展矩阵的列向量进行并行处理，其中，对于所述扩展矩阵的每一列列向量上的数据采用串行计算的方式进行计算，得到计算结果，所述计算结果包括所述目标矩阵的逆矩阵和与所述逆矩阵具有相同大小的单位矩阵；以及

输出单元，用于将所述计算结果输出至所述中央处理器，其中，所述中央处理器从所述计算结果中提取所述逆矩阵。

7.根据权利要求6所述的矩阵求逆处理装置，其特征在于，所述第一建立单元包括：

确定模块，用于根据所述扩展矩阵的大小确定分割基数；

划分模块，用于按照所述分割基数对所述扩展矩阵的行向量和列向量进行划分，得到多个数据段；以及

构建模块，用于根据所述分割基数和所述多个数据段的数量构建线程块结构，形成所述全局网格结构，其中，所述全局网格结构包括与所述多个数据段一一对应的线程块，每个线程块具有与所述分割基数相同数量的线程。

8.根据权利要求6所述的矩阵求逆处理装置，其特征在于，所述计算单元包括：

第一计算模块，用于计算所述扩展矩阵中当前行向量的系数向量，所述系数向量包括所述扩展矩阵中所述当前行向量之外的其他行向量对于所述当前行向量的系数；

第二计算模块，用于利用所述全局网格结构中线程计算该线程在所述扩展矩阵的映射位置上的变换结果，用所述变换结果替换在所述扩展矩阵的映射位置上的数据，得到替换后的扩展矩阵；

判断模块，用于判断所述当前行向量是否是所述扩展矩阵的最后一行行向量；

所述第一计算模块还用于如果判断出所述当前行向量不是所述扩展矩阵的最后一行行向量，则将所述当前行向量的下一行行向量作为所述当前行向量，计算当前行向量的系数向量；以及

处理模块，用于如果判断出所述当前行向量是所述扩展矩阵的最后一行行向量，则对所述替换后的扩展矩阵进行单位化处理，得到所述计算结果。

9.根据权利要求8所述的矩阵求逆处理装置，其特征在于，所述第一计算模块包括：

第一获取子模块，用于获取所述当前行向量的对角位置数据；

第二获取子模块，用于获取所述对角位置数据所在的列向量上的数据；

除法子模块，用于依次将对角位置数据与所述对角位置数据所在的列向量上的数据相除，得到所述系数向量。

10.根据权利要求8所述的矩阵求逆处理装置，其特征在于，所述矩阵求逆处理装置还包括：

第二建立单元，用于在利用所述全局网格结构中线程计算该线程在所述扩展矩阵的映射位置上的变换结果之前，建立共享显存空间，所述共享显存空间用于存储所述当前行向量的数据；

所述第二计算模块还用于所述全局网格结构中线程从所述共享显存空间调用所述当前行向量的数据用于计算该线程在所述扩展矩阵的映射位置上的变换结果。