CN102880594B - 基于多核dsp的并行矩阵全选主元高斯约旦求逆方法 - Google Patents

Abstract

本发明提供了一种基于多核DSP的并行矩阵全选主元高斯约旦求逆算法，解决了现阶段关于求逆矩阵的算法不能满足所有矩阵，并且现有的矩阵求逆算法已经不能满足高速运算要求的缺陷。本发明将矩阵A划分为p块子矩阵后依此发送给1～（p-1）号处理器，拥有当前主行元素的处理器对主行元素进行更新，并发送给除第0号处理器和自身处理器外其余的处理器后，1～（p-1）号处理器一样对自身的子矩阵进行更新，并分别将自身的子矩阵回传给第0号处理器，然后选择上一主行元素的下一行元素作为新的主行元素重复上述过程，直到主行元素遍历完原矩阵的所有对角线元素。通过上述方案，本发明达到了既能针对所有矩阵可实现求逆功能、又满足高速运算要求的效果。

Description

基于多核DSP的并行矩阵全选主元高斯约旦求逆方法

技术领域

本发明属于移动通信领域，具体的讲，涉及一种基于多核DSP的并行矩阵全选主元高斯约旦求逆方法。

背景技术

矩阵求逆是工程实践中常用也是较为繁琐的运算过程，现有的普通矩阵求你算法种类繁多，其中绝对大部分的算法复杂程度可以满足一般的工程应用要求。然而随着应用的日益复杂，这些的矩阵求逆算法已经不能满足高速运算的要求。

因此现有技术中为了满足矩阵求逆高速运算的要求，提出了一些高速的求逆算法，但是这些算法都是针对一些特殊矩阵而实现，如针对三对角线矩阵、五对角线矩阵和三角矩阵逆矩阵的算法，这些算法对矩阵的类型有特定的要求，因此开发出既能针对所有矩阵均可实现求逆功能、并且满足高速运算要求的求逆算法已成为现阶段人们研究的重点。

发明内容

本发明的目的在于克服现有技术中关于求逆矩阵的算法不能满足所有矩阵，并且现有的矩阵求逆算法已经不能满足高速运算要求的缺陷，提供一种基于多核DSP的并行矩阵全选主元高斯约旦求逆方法，充分地发挥多核DSP的高性能优势，解决上述缺陷。

本发明采用的技术方案如下：

基于多核DSP的并行矩阵全选主元高斯约旦求逆方法，通过以下步骤实现：

(1)共有p块处理器，将处理器依次编号为0～(p-1)，并通过第0号处理器对原矩阵A进行全选主元，得到主行元素，其中p为正整数；

(2)将进行了全选主元后的矩阵A划分为p块子矩阵，第0号处理器保留其中一个子矩阵，并将另外p-1块子矩阵分别发送给第1～(p-1)号处理器，使每个处理器分别接收一块子矩阵；

(3)拥有当前主行元素的处理器对主行元素进行更新，并将更新后的主行元素发送给除第0号处理器和自身处理器外其余的处理器；

(4)除第0号处理器和自身处理器外其余的处理器收到更新后的主行元素数值后，第1～(p-1)号处理器分别对自身的子矩阵进行更新；之后分别将自身的子矩阵回传给第0号处理器；同时，子矩阵中拥有当前主行的处理器也利用相同的公式对其子矩阵中非主行元素进行更新；步骤(4)中1～(p-1)号处理器使用公式：

${\mathrm{\α}}_{\mathrm{ij}}={\mathrm{\α}}_{\mathrm{ij}}^{*}-{\mathrm{\α}}_{\mathrm{ik}}^{*}\×{\mathrm{\α}}_{\mathrm{kj}}$

和公式：

${\mathrm{\α}}_{\mathrm{ik}}=-{\mathrm{\α}}_{\mathrm{ik}}^{*}\×{\mathrm{\α}}_{\mathrm{kk}}$

对自身的子矩阵进行更新，其中i为整数且0≤i＜n，同时i≠k；

(5)第0号处理器进行整理，再进行全选主元；

(6)选择原矩阵A主行元素的下一行元素作为新的主行元素，重复步骤(3)～(5)，直到主行元素遍历完原矩阵的所有对角线元素，此时第0号处理器整理得到的矩阵即为原矩阵A的逆矩阵。

为了便于查询，所述步骤(1)和步骤(6)中对原矩阵A进行全选主元后，记录交换元素的行列号。

为了分散计算，所述步骤(2)中具体实现步骤如下：

(2a)从原矩阵中的第0行开始，依次将原矩阵的第j行分配给第(jmod

p)号的处理器，其中所述j为整数且j≥0；

(2b)每一个处理器中被分配的所有行组成一个子矩阵，在所有子矩阵形成后，第0号处理器保留其中一个子矩阵，并将另外p-1块子矩阵分别发送给第1～(p-1)号处理器，使每个处理器分别接收一块子矩阵。

再进一步的，所述步骤(3)中子矩阵中拥有当前主行的处理器利用公式：

${\mathrm{\α}}_{\mathrm{kk}}=1/{\mathrm{\α}}_{\mathrm{kk}}^{*}$

和公式：

${\mathrm{\α}}_{\mathrm{kj}}={\mathrm{\α}}_{\mathrm{kj}}^{*}\×{\mathrm{\α}}_{\mathrm{kk}}$

对主行元素进行更新，其中式中α*为子矩阵中元素值，α为更新后子矩阵中的元素值，n为子矩阵行数，j为整数且0≤j＜n，k为整数且0≤k＜n，同时j≠k。

为了满足不同数据传输要求，所述步骤(3)中拥有当前主行的处理器利用公式对主行元素进行更新后，使用EDMA3或SRIO将更新后的主行元素发送给第1～(p-1)号处理器中除自身处理器外其余的处理器。

为了再次满足不同数据传输要求，所述步骤(5)中第1～(p-1)号处理器中的子矩阵完成更新后，使用EDMA3或SRIO将其子矩阵回传给第0号处理器。

与现有技术相比，本发明具有以下有益效果：

(1)本发明利用多核DSP平台进行矩阵的求逆计算，即将矩阵划分为多个子矩阵并将每个子矩阵放入每一个处理器，每一个处理器对于自身子矩阵再分别进行更新，这样分散了矩阵求逆计算，可以有效地提高矩阵求逆运算的速度；通过实验证明，在矩阵阶数越大，并且并行运算的处理器个数越多时，运算速度越快，效果越明显。

(2)本发明通过全选主元高斯约旦算法达到了占用内存小的同时求解精度高、运算速度快的效果，高斯约旦算法主要为求解矩阵的逆矩阵，采用多核DSP平台通过高斯约旦算法求解矩阵的逆矩阵，结合了多核平台与高斯约旦算法共同的优点，达到了事半功倍的效果。因为其不仅大大提高了求解的速度，同时还提高了矩阵的求解精度，质量和速度同时得到了保证。

(3)本发明采用的是多核DSP，DSP是一种独特的微处理器，是以数字信号来处理大量信息的器件。它不仅具有可编程性，而且实时运行速度可达每秒数以千万条复杂指令程序，远远超过通用微处理器，是数字化电子世界中日益重要的电脑芯片。它的强大数据处理能力和高运行速度，是最值得称道的两大特色，多核DSP被应用于本发明，在很大程度上提高了相互通信的处理器之间的数据交换速度。

(4)本发明通过行交叉划分将矩阵A划分为p块子矩阵，即将原矩阵的第j行分配给第(jmodp)号的处理器，再由分配该处理器所有原矩阵的行组成一个子矩阵，在分配好所有的子矩阵后，将其中的p-1块子矩阵发送给对应的第1～(p-1)号处理器。该方法被称为行交叉法，通过上述方法来划分矩阵A的p块子矩阵，不仅简便快捷，充分发挥了多核处理器的优势，同时还保证了后面步骤的正常执行。

(5)本发明还使用了EDMA3或SRIO将子矩阵从第0号处理器发送给其余处理器，或者将其余处理器的子矩阵回传给第0号处理器，EDMA3或SRIO能满足不同的数据传输需求，是未来十几年中嵌入式系统互联的最佳选择之一。

附图说明

图1为本发明—实施例1中并行算法流程图。

图2为本发明—实施例1中主行元素位于第0号处理器子矩阵时的通信示意图。

图3为本发明—实施例2中主行元素位于第1号处理器子矩阵时的通信示意图。

图4为本发明—实施例3中主行元素位于第2号处理器子矩阵时的通信示意图。

图5为本发明—实施例4中主行元素位于第3号处理器子矩阵时的通信示意图。

具体实施方式

下面结合附图与实施例对本发明作进一步说明，本发明的实施方式包括但不限于下列实施例。

实施例1

如图1所示的基于多核DSP的并行矩阵求逆算法，该算法为全选主元高斯约旦算法，占用内存小，运算速度快且精度较高，易于并行实现等。

高斯约旦求逆算法矩阵元素的每次更新都是以对角线上的元素为基础，分别对该元素同行不同列，同列不同行以及既不同列也不同行的元素以不同的公式进行更新。其中每一次更新中当前对角线元素所在的行称为主行，主行上的各矩阵元素称为主行元素，其余矩阵元素称为非主行元素。

如图1、2所示，本实施例中使用j＝5块处理器对n阶矩阵A求逆，处理器从0到4进行编号。其中数字信号处理器DSP1中包括第0号处理器、第1号处理器和第2号处理器，为三核处理器；数字信号处理器DSP2包括第3号处理器和第4号处理器，为双核处理器。图2为本实施例中主行元素位于第0号处理器子矩阵时的通信示意图。具体的并行算法如下：

第1步：第0号处理器对原矩阵A进行全选主元，并记录交换元素的行列号。

其中所述全选主元即主列消元，并得到主行元素，这里的主行为矩阵A的第一行。

第2步：将矩阵A划分为5个子矩阵，其中具体实现步骤如下：

(2a)从矩阵A中的第一行开始，依次将原矩阵的每一行进行分配，即将第j行分配给第(jmod5)号的处理器，其中所述j为正整数；

其中mod为取余运算，即计算j除以5的余数；

(2b)每一个处理器中被分配的所有行组成一个子矩阵，在所有子矩阵形成后，将其中的4块子矩阵发送给对应的第1～4号处理器。

第3步：其中子矩阵中拥有主行的处理器为第0号处理器，因此第0号处理器利用公式：

${\mathrm{\α}}_{\mathrm{kk}}=1/{\mathrm{\α}}_{\mathrm{kk}}^{*}$

和公式：

${\mathrm{\α}}_{\mathrm{kj}}={\mathrm{\α}}_{\mathrm{kj}}^{*}\×{\mathrm{\α}}_{\mathrm{kk}}$

对主行元素进行更新，并使用EDMA3将更新后的主行元素发送给第1号处理器、第2号处理器，使用SRIO将更新后的主行元素发送给第3号处理器、第4号处理器，其中式中α*为子矩阵中元素值，α为更新后子矩阵中的元素值，n为子矩阵行数，j为整数且0≤j＜n，k为整数且0≤k＜n，同时j≠k。

第4步：其余处理器根据接收到更新后的主行元素数值后，使用公式

${\mathrm{\α}}_{\mathrm{ij}}={\mathrm{\α}}_{\mathrm{ij}}^{*}-{\mathrm{\α}}_{\mathrm{ik}}^{*}\×{\mathrm{\α}}_{\mathrm{kj}}$

和公式

${\mathrm{\α}}_{\mathrm{ik}}=-{\mathrm{\α}}_{\mathrm{ik}}^{*}\×{\mathrm{\α}}_{\mathrm{kk}}$

对其子矩阵进行更新，其中i为整数且0≤i＜n，同时i≠k。

第5步：第1～4号处理器的子矩阵完成更新后，第1号处理器、第2号处理器使用EDMA3，同时第3号处理器、第4号处理器使用SRIO将其子矩阵回传给第0号处理器。

第6步：第0号处理器对来自其他处理器的矩阵进行整理，得到一个新的完整的矩阵，再对其进行全选主元，并记录交换元素的行列号。

第7步：选择主行的下一行作为新的主行，重复第三到第六步，直到主元遍历完原矩阵的所有对角线元素，此时第0号处理器整理得到矩阵就是原矩阵的逆矩阵。

实施例2

与实施例1的区别在于步骤3～5不同。

本实施例步骤3～5具体实现方式如下：

第3步：其中子矩阵中拥有主行的处理器第1号处理器，因此第1号处理器利用公式：

${\mathrm{\α}}_{\mathrm{kk}}=1/{\mathrm{\α}}_{\mathrm{kk}}^{*}$

和公式：

${\mathrm{\α}}_{\mathrm{kj}}={\mathrm{\α}}_{\mathrm{kj}}^{*}\×{\mathrm{\α}}_{\mathrm{kk}}$

对主行元素进行更新，并使用EDMA3将更新后的主行元素发送给第2号处理器，使用SRIO将更新后的主行元素发送给第3号处理器、第4号处理器，其中式中α*为子矩阵中元素值，α为更新后子矩阵中的元素值，n为子矩阵行数，j为整数且0≤j＜n，k为整数且0≤k＜n，同时j≠k，如图3所示，该图为主行元素位于第1号处理器子矩阵时的通信示意图。

第4步：其余处理器根据接收到更新后的主行元素数值，使用公式：

${\mathrm{\α}}_{\mathrm{ij}}={\mathrm{\α}}_{\mathrm{ij}}^{*}-{\mathrm{\α}}_{\mathrm{ik}}^{*}\×{\mathrm{\α}}_{\mathrm{kj}}$

和公式：

${\mathrm{\α}}_{\mathrm{ik}}=-{\mathrm{\α}}_{\mathrm{ik}}^{*}\×{\mathrm{\α}}_{\mathrm{kk}}$

对其子矩阵进行更新，其中i为整数且0≤i＜n，同时i≠k。

第5步：第2～4号处理器的子矩阵完成更新后，第1号处理器、第2号处理器使用EDMA3，同时第3号处理器、第4号处理器使用SRIO将其子矩阵回传给第0号处理器。

其余步骤与实施例1相同。

实施例3

与实施例1的区别在于步骤3～5不同。

本实施例步骤3～5具体实现方式如下：

第3步：其中子矩阵中拥有主行的处理器为第2号处理器，因此第2号处理器利用公式：

${\mathrm{\α}}_{\mathrm{kk}}=1/{\mathrm{\α}}_{\mathrm{kk}}^{*}$

和公式：

${\mathrm{\α}}_{\mathrm{kj}}={\mathrm{\α}}_{\mathrm{kj}}^{*}\×{\mathrm{\α}}_{\mathrm{kk}}$

对主行元素进行更新，并使用EDMA3将更新后的主行元素发送给第1号处理器，使用SRIO将更新后的主行元素发送给第3号处理器、第4号处理器。

第4步：其余处理器根据接收到更新后的主行元素数值后，使用公式

${\mathrm{\α}}_{\mathrm{ij}}={\mathrm{\α}}_{\mathrm{ij}}^{*}-{\mathrm{\α}}_{\mathrm{ik}}^{*}\×{\mathrm{\α}}_{\mathrm{kj}}$

和公式

${\mathrm{\α}}_{\mathrm{ik}}=-{\mathrm{\α}}_{\mathrm{ik}}^{*}\×{\mathrm{\α}}_{\mathrm{kk}}$

对其子矩阵进行更新，其中i为整数且0≤i＜n，同时i≠k。

第5步：第1号处理器、第3号处理器、第4号处理器的子矩阵完成更新后，第1号处理器、第2号处理器使用EDMA3，同时第3号处理器、第4号处理器使用SRIO将其子矩阵回传给第0号处理器，其中式中α*为子矩阵中元素值，α为更新后子矩阵中的元素值，n为子矩阵行数，j为整数且0≤j＜n，k为整数且0≤k＜n，同时j≠k。如图4所示，该图为主行元素位于第2号处理器子矩阵时的通信示意图。

其余步骤与实施例1相同。

实施例4

与实施例1的区别在于步骤3～5不同。

本实施例步骤3～5具体实现方式如下：

第3步：其中子矩阵中拥有主行的处理器为第3号处理器，因此第3号处理器利用公式：

${\mathrm{\α}}_{\mathrm{kk}}=1/{\mathrm{\α}}_{\mathrm{kk}}^{*}$

和公式：

${\mathrm{\α}}_{\mathrm{kj}}={\mathrm{\α}}_{\mathrm{kj}}^{*}\×{\mathrm{\α}}_{\mathrm{kk}}$

对主行元素进行更新，并使用SRIO将更新后的主行元素发送给第1号处理器、第2号处理器，使用EDMA3将更新后的主行元素发送给第4号处理器。

第4步：其余处理器根据接收到更新后的主行元素数值后，使用公式

${\mathrm{\α}}_{\mathrm{ij}}={\mathrm{\α}}_{\mathrm{ij}}^{*}-{\mathrm{\α}}_{\mathrm{ik}}^{*}\×{\mathrm{\α}}_{\mathrm{kj}}$

和公式

${\mathrm{\α}}_{\mathrm{ik}}=-{\mathrm{\α}}_{\mathrm{ik}}^{*}\×{\mathrm{\α}}_{\mathrm{kk}}$

对其子矩阵进行更新，其中i为整数且0≤i＜n，同时i≠k。

第5步：从第1～4号处理器的子矩阵完成更新后，第1号处理器、第2号处理器使用SRIO同时第4号处理器使用EDMA3将其子矩阵回传给第0号处理器，其中式中α*为子矩阵中元素值，α为更新后子矩阵中的元素值，n为子矩阵行数，j为整数且0≤j＜n，k为整数且0≤k＜n，同时j≠k。如图5所示，该图为主行元素位于第3号处理器子矩阵时的通信示意图。

根据本发明的内容，将并行的高斯约旦求逆算法在TMS320C6474ENM板上进行了实现，得出的逆矩阵结果正确。

表1高斯约旦算法在不同平台下占用的CPU时钟数

矩阵阶数	12	24	36	48
					单核	1,224,192	9,454,188	31,486,992	74,354,646
两核	1,242,591	9,640,287	32,343,931	77,696,788
					四核	1,098,168	8,089,545	26,910,182	66,073,290

其中，TMS320C6474DSP的CPU时钟频率为1GHz。

通过上述表格得到，当参与运算的处理器个数较大以及矩阵的阶数较大时，算法的效果非常好。

按照上述实施例，便可很好的实现本发明。

Claims (6)

1.基于多核DSP的并行矩阵全选主元高斯约旦求逆方法，其特征在于，通过以下步骤实现：

(1)共有p块处理器，将处理器依次编号为0～(p-1)，并通过第0号处理器对原矩阵A进行全选主元，得到主行元素，其中p为正整数；

(2)将进行了全选主元后的矩阵A划分为p块子矩阵，第0号处理器保留其中一个子矩阵，并将另外p-1块子矩阵分别发送给第1～(p-1)号处理器，使每个处理器分别接收一块子矩阵；

(3)拥有当前主行元素的处理器对主行元素进行更新，并将更新后的主行元素发送给除第0号处理器和自身处理器外其余的处理器；

(4)除第0号处理器和自身处理器外其余的处理器收到更新后的主行元素数值后，第1～(p-1)号处理器分别对自身的子矩阵进行更新；之后分别将自身的子矩阵回传给第0号处理器，第1～(p-1)号处理器使用公式：

${\mathrm{\α}}_{ij}={\mathrm{\α}}_{ij}^{*}-{\mathrm{\α}}_{ik}^{*}\×{\mathrm{\α}}_{kj}$

和公式：

${\mathrm{\α}}_{ik}=-{\mathrm{\α}}_{ik}^{*}\×{\mathrm{\α}}_{kk}$

对自身的子矩阵进行更新，其中i为整数且0≤i＜n，j为整数且0≤j＜n，k为整数且0≤k＜n，同时i≠k，α*为子矩阵中元素值，α为更新后子矩阵中的元素值，n为子矩阵行数；同时，子矩阵中拥有当前主行的处理器也利用相同的公式对其子矩阵中非主行元素进行更新；

(5)第0号处理器进行整理，再进行全选主元；

(6)选择原矩阵A主行元素的下一行元素作为新的主行元素，重复步骤(3)～(5)，直到主行元素遍历完原矩阵的所有对角线元素，此时第0号处理器整理得到的矩阵即为原矩阵A的逆矩阵。

2.根据权利要求1所述的基于多核DSP的并行矩阵全选主元高斯约旦求逆方法，其特征在于，所述步骤(1)和步骤(6)中对原矩阵A进行全选主元后，记录交换元素的行列号。

3.根据权利要求2所述的基于多核DSP的并行矩阵全选主元高斯约旦求逆方法，其特征在于，所述步骤(2)中具体实现步骤如下：

(2a)从原矩阵中的第0行开始，依次将原矩阵的第j行分配给第(jmodp)号的处理器，其中所述j为整数且j≥0；

(2b)每一个处理器中被分配的所有行组成一个子矩阵，在所有子矩阵形成后，第0号处理器保留其中一个子矩阵，并将另外p-1块子矩阵分别发送给第1～(p-1)号处理器，使每个处理器分别接收一块子矩阵。

4.根据权利要求3所述的基于多核DSP的并行矩阵全选主元高斯约旦求逆方法，其特征在于，所述步骤(3)中子矩阵中拥有当前主行的处理器利用公式：

${\mathrm{\α}}_{kk}=1/{\mathrm{\α}}_{kk}^{*}$

和公式：

${\mathrm{\α}}_{kj}={\mathrm{\α}}_{kj}^{*}\×{\mathrm{\α}}_{kk}$

对主行元素进行更新，其中式中α*为子矩阵中元素值，α为更新后子矩阵中的元素值，n为子矩阵行数，j为整数且0≤j＜n，k为整数且0≤k＜n，同时j≠k。

5.根据权利要求4所述的基于多核DSP的并行矩阵全选主元高斯约旦求逆方法，其特征在于，所述步骤(3)中拥有当前主行的处理器利用公式对主行元素进行更新后，使用EDMA3或SRIO将更新后的主行元素发送给第1～(p-1)号处理器中除自身处理器外其余的处理器。

6.根据权利要求1～5任意一项所述的基于多核DSP的并行矩阵全选主元高斯约旦求逆方法，其特征在于，所述步骤(5)中第1～(p-1)号处理器中的子矩阵完成更新后，使用EDMA3或SRIO将其子矩阵回传给第0号处理器。