CN102662917B - 正定Hermite矩阵Cholesky分解高速脉动阵列的设计方法 - Google Patents

Abstract

本发明公开了正定Hermite矩阵的Cholesky分解高速脉动阵列的设计方法。本发明将处理单元划分为圆形处理单元和方形处理单元，并对不同类别的处理单元进行了功能划分，保证硬件设计具有规律性，实现数据的并行化流水处理。所述设计方法适合所有维数为M×M的正定Hermite矩阵的Cholesky分解，具有通用性，适用范围广。

Description

正定Hermite矩阵Cholesky分解高速脉动阵列的设计方法

技术领域

本发明涉及正定Hermite矩阵的分解方法，具体是指正定Hermite矩阵Cholesky分解脉动阵列的方法。

背景技术

在科学计算和工程领域中，Cholesky分解是一个非常重要的线性代数法。尤其是在第三代移动通信TD-SCDMA系统中，多种联合检测算法都要涉及到系统矩阵A求逆问题，而在对矩阵A求逆时，首先要对矩阵A进行分解，当用户比较多且接收天线比较多时，将矩阵A分解的计算量将会比较庞大，复杂度也比较高。Cholesky分解因为利用了Hermite矩阵的性质，在对矩阵A进行分解时，有效地降低了运算复杂度。

而基于Cholesky分解正定Hermite矩阵的方法大多是根据Cholesky分解公式直接编写的，这种方法因为把分解当作一个整体，没有体现分块设计思想，特别是矩阵维数较大时，数据之间的关系复杂而混乱，难以保证设计的正确性，并且因为数据计算较复杂，因此硬件的并行处理速度并不高。

发明内容

本发明的目的在于克服通过Cholesky分解公式直接编写的硬件设计不适用于维数较大的矩阵、并且硬件并行处理速度不高的缺陷，提供一种涉及M×M维的正定Hermite矩阵的Cholesky分解高速脉动阵列的设计方法。

本发明的目的通过下述技术方案实现：

正定Hermite矩阵Cholesky分解高速脉动阵列的设计方法，包括以下步骤：

(a)根据维数为M×M正定Hermite矩阵确定脉动阵列的维数以及构架；

(b)根据脉动阵列的架构对脉动阵列的处理单元进行分类；

(c)根据处理单元类别设计处理单元结构；

(d)根据脉动阵列的架构连接处理单元；

(e)输出脉动阵列。同时，所述步骤(a)中脉动阵列为M×M维的包含斜边的上三角处理单元阵列，且所述脉动阵列共包含(1+M)M/2个处理单元。

并且，所述步骤(a)中脉动阵列为M×M维的上三角单元处理阵列，且所述脉动阵列共有(1+M)M/2个处理单元。

进一步的，所述步骤(c)中对圆形处理单元结构的设计具体包括以下步骤：

(c1)将脉动阵列上位于第M列、第M行的处理单元划分为第一类圆形处理单元，所述第一类圆形处理单元内置有一个开方器和一个寄存器，且包含一个输入端和一个输出端；

(c2)将脉动阵列斜边上的其他处理单元划分为第二类圆形处理单元，所述第二类圆形处理单元内置有一个开方器、一个除法器和两个寄存器，且包含一个输入端和两个输出端。

更进一步的，所述步骤(c)中对方形处理单元结构的设计包括以下步骤：

(I)将平行于圆形处理单元且与圆形处理单元相邻的处理单元划分为第一类方形处理单元，所述第一类方形处理单元由两个乘法器、一个加法器、一个选择器、两个移位寄存器和两个寄存器组成，且包含两个输入端和两个输出端；

(II)将平行于第一类方形处理单元且与第一类方形处理单元相邻的处理单元划分为第二类方形处理单元，所述第二类方形处理单元由两个乘法器、一个加法器、三个选择器两个移位寄存器和两个寄存器组成，且包含三个输入端和两个输出端；

(III)以此类推，将平行于第K-1类方形处理单元且与第K-1类方形处理单元相邻的处理单元划分为第K类方形处理单元，所述第K类方形处理单元由两个乘法器、一个加法器、三个选择器、两个移位寄存器和两个寄存器组成，其中第K类方形处理单元有K+1个输入端，两个输出端。

并且所述K≤M。

为了更好的实现本发明，所述步骤(d)中在连接处理单元时，脉动阵列上同时连接有延迟单元。

本发明较现有技术相比，具有以下优点及有益效果：

(1)本发明采用脉动阵列结构，数据按预先确定的方式在脉动阵列的处理单元间有规律地传输，所有处理单元能同时并行地对流经它的数据进行处理，解决了常规设计中硬件并行处理速度不高的缺陷，提高了处理单元数据处理效率；

(2)本发明对整个矩阵进行了分解处理，将脉动阵列结构划分为不同类型的处理单元，实现处理单元的分工合作，保证硬件设计具有规律性，不仅适合维数较小的矩阵，同时适合维数较大的矩阵计算，从而在提高计算正确率的同时，扩大了适用范围，使其具有了更为广阔的应用领域；

(3)本发明在方形处理单元中尽量复用乘法器，避免了方形处理单元使用乘法器数量随分解矩阵维数的增加而快速增加的缺陷，因而有效的减少了硬件资源的开销；

(4)本发明中，将数据输入处理单元阵列后，处理单元阵列能自动完成数据处理，脉动阵列的运行速度大大提高，同时简化了逻辑控制；

(5)本发明的处理单元结构和阵列简单，并且具有一致的规则，可以达到很高的模块化程度，非常适合超大规模集成电路的设计和制造；

(6)本发明处理单元的规则性强，适合所有维数为M×M的正定Hermite矩阵的Cholesky分解，具有通用性，适用范围广。

附图说明

图1为本发明—实施例4×4维的脉动阵列结构图。

图2为本发明—实施例的输入输出时序图。

图3为本发明的第一类圆形处理单元的结构图。

图4为本发明的第二类圆形处理单元的结构图。

图5为本发明的第一类方形处理单元的结构图。

图6为本发明的第二类方形处理单元的结构图。

图7为本发明的第三类方形处理单元的结构图。

图8为本发明—实施例的Verilog HDL定点仿真与Matlab浮点仿真的误差曲线图。

图9为本发明—实施例的流程图。

具体实施方式

下面结合实施例对本发明作进一步的详细说明，但本发明的实施方式不限于此。

实施例

本实施例中，使用Matlab仿真Cholesky分解脉动阵列实现正定Hermite矩阵的分解。

正定Hermite矩阵A为4×4维矩阵，由Hermite矩阵的性质：A＝A^H，即矩阵的元素将A表示为：

$A=\left(\begin{array}{cccc}{a}_{11}& a_{21}^{*}& a_{31}^{*}& a_{41}^{*}\\ a_{21}& a_{22}& a_{32}^{*}& a_{42}^{*}\\ a_{31}& a_{32}& a_{33}& a_{43}^{*}\\ a_{41}& a_{42}& a_{43}& a_{44}\end{array}\right)$

显然,矩阵A满足Cholesky分解的要求，因此待分解矩阵处于斜边上半部分的元素完全可以用斜边左下方的元素表示，输入矩阵A包含斜边元素在内的下三角的元素就足以表示整个矩阵A。

下三角矩阵L中每个元素l_ij(i＝1,2,3,4；j＝1,2,3,4)的计算公式如下：

$l_{\mathrm{ij}}=\left\{\begin{array}{ccc}\sqrt{a_{\mathrm{ii}}-\underset{k=1}{\overset{i=1}{\mathrm{\&Sigma;}}}{l}_{\mathrm{ik}}{l}_{\mathrm{ik}}^{*}},& \mathrm{for}& i=j\\ \frac{1}{l_{\mathrm{jj}}}(a_{\mathrm{ij}}-\underset{k=1}{\overset{j-1}{\mathrm{\&Sigma;}}}{l}_{\mathrm{ik}}{l}_{\mathrm{jk}}^{*})& \mathrm{for}& i>j\\ 0,& \mathrm{for}& i<j\end{array}\right.$

其中，下三角矩阵L的主对角元素l_ij均为正实数。

根据上述公式的展开可得矩阵L第一列的所有元素：

$\left\{\begin{array}{c}{l}_{11}=\sqrt{a_{11}}\\ l_{21}=\frac{a_{21}}{l_{11}}\\ l_{31}=\frac{a_{31}}{l_{11}}\\ l_{41}=\frac{a_{41}}{l_{11}}\end{array}\right.$

分析上式可知，作为一个模块，包含开方和除法两个部分，而该模块的输出l₁₁作为l₂₁、l₃₁、l₄₁的输入，因此便可以得到矩阵L第一列的所有元素。

矩阵L第二列的计算公式如下：

$\left\{\begin{array}{c}{l}_{22}=\sqrt{a_{22}-l_{21}{l}_{21}^{*}}\\ l_{32}=\frac{1}{l_{22}}(a_{32}-l_{31}{l}_{21}^{*})\\ l_{42}=\frac{1}{l_{22}}(a_{42}-l_{41}{l}_{21}^{*})\end{array}\right.$

矩阵L第三列的计算公式如下：

$\left\{\begin{array}{c}{l}_{33}=\sqrt{a_{33}-l_{31}{l}_{31}^{*}-l_{32}{l}_{32}^{*}}\\ l_{43}=\frac{1}{l_{43}}(a_{43}-l_{41}{l}_{31}^{*}-l_{42}{l}_{32}^{*})\end{array}\right.$

根据第二列和第三列的计算公式可得，也同样包含开方和除法两部分。矩阵L第一列的计算结果l₂₁、l₃₁、l₄₁都作为第二列、第三列计算的输入，具体说来，l₂₁作为求l₂₂的输入，l₃₁作为求l₃₂的输入，l₄₁作为求l₄₂的输入。

同时，第四行第四列的元素计算公式如下：

$l_{44}=\sqrt{a_{44}-l_{41}{l}_{41}^{*}-l_{42}{l}_{42}^{*}-l_{43}{l}_{43}^{*}}$

基于以上数据的流动规则，可得上述元素计算公式的脉动阵列结构图，如图1所示，本发明—实施例的脉动阵列是一个4×4维的包含斜边在内的上三角处理单元阵列，所述脉动阵列共包含10个处理单元，处理单元简称PE。

根据图1所示，脉动阵列上的PE分为圆形PE和方形PE，根据上述公式规律可得，脉动阵列斜边上的PE为圆形PE，其余为方形PE。

本发明—实施例旨在实现4×4维正定Hermite矩阵通过Cholesky分解高速脉动阵列的方法实现矩阵的分解，输入矩阵A经过一定处理后按每一行并行输出其分解矩阵L。

对应的输入输出时序如图2所示，在输出时序图2后，如图9所示，由下列步骤实现本发明—实施例：

步骤一：确定脉动阵列的维数以及构架

其中矩阵H是均值为0、方差为1的4阶随机复矩阵，由H生成4×4维正定Hermite矩阵A，通过Cholesky分解设计出其脉动阵列，并用FPGA实现，写成MATLAB代码为：

n＝4；

H＝(randn(n,n)+1j*randn(n,n))/sqrt(2)；

A＝H’*H；

用(y1，y2，y3)表示一个定点数，其中y1为符号位，取值为0时表示无符号，取值为1时表示有符号，y2表示整数部分的位宽，y3表示小数部分的位宽。在本发明—实施例中，将输入矩阵A定点为(1，5，19)，即表示矩阵为有符号位，5位整数位和19位小数位。

根据矩阵L每个元素的计算公式确定脉动阵列的PE构架，并根据矩阵元素的表示方式对各个PE进行编号，如图1所示。

步骤二：对圆形PE结构及功能进行设计。

斜边上PE为圆形PE。其中PE₄₄只包含开方运算，因此被划分为第一类圆形PE；斜边上其余处理单元包含开方以及除法运算，因此被划分为第二类圆形处理单元，包括PE₁₁、PE₂₂、PE₃₃。

第一类圆形处理单元结构如图3所示，所述第一类圆形PE有一个输入、一个输出。该类圆形PE输入数据rx_yin，然后进行开方，开方的结果作为该圆形PE的第一个输出tx_xout，计算公式如下：

$\mathrm{tx}\_\mathrm{xout}=\sqrt{\mathrm{rx}\_\mathrm{yin}}$

第二类圆形PE结构如图4所示，所述第二类圆形PE有一个输入和两个输出，输入数据rx_yin进行开方后输出结果tx_xout，同时tx_xout进行求倒数，计算公式如下：

$\mathrm{tx}\_\mathrm{xout}\_d=\frac{1}{\mathrm{tx}\_\mathrm{xout}}$

步骤三：对方形PE的结构和功能进行设计。

方形PE被划分为：第一类方形PE，包括PE₁₂、PE₂₃、PE₃₄；第二类方形PE，包括PE₁₃、PE₂₄；第三类方形PE，包括PE₁₄。

第一类方形PE的结构框图如图5所示，该类方形PE的输入rx_xin连接的是与它处于同一行的第一类圆形PE的输出tx_xout_d；所述第一类圆形PE一定位于第一类方形PE的左侧，同时所述第一类方形PE的输入rx_yin与它处于同一列的第二类方形PE的输出tx_xout相连。输入端rx_xin每M个时钟输入一个数据；输入端rx_yin每M个时钟输入两个数据，第一个数据为rx_yin(1)，第二个数据为rx_yin(2)。rx_xin与rx_yin(1)作乘法，然后将乘积作为tx_xout输出，计算公式如下：

tx_xout＝rx_yin(1)×rx_xin

输出的同时，将tx_xout反馈回另一个乘法器的输入端，同时另一个输入端输入tx_xout的共轭tx_xout*，然后将该乘法器输出的结果取反后与rx_yin(2)作加法运算，计算公式如下：

tx_yout＝-tx_xout^*×tx_xout+rx_yin(2)

因此第一类方形PE由两个乘法器、一个加法器和一个选择器组成，计算过程中复用了一次乘法器，减少硬件资源的开销。

第二类方形PE的结构框图如图6所示，该类方形PE的输入rx_xin1连接的是与它处于同一行的第二类圆形PE的输出tx_xout_d，输入rx_xin2连接的是与它处于同一行的第一个方形PE的输出tx_xout，另一个输入rx_yin与它处于同一列上方的第三类方形PE的输出tx_yout相连接，其中输入端rx_xin1与输入端rx_xin2每M个时钟输入一个数据，输入端rx_yin每M个时钟输入三个数据，分别为rx_yin(1)，rx_yin(2)，rx_yin(3)。

rx_xin1的第一个数据与rx_yin的第一个数据rx_yin(1)作乘法，然后将乘积作为tx_xout输出，计算公式如下：

tx_xout＝rx_yin(1)×rx_xin1

输出的同时，将tx_xout反馈回另一个选择器的输入端，一个选择器的输入端输入全部tx_xout,另一个选择器的一个输入端输入tx_xout*，选择器将rx_xin2的共轭rx_xin2*以及tx_xout输入到另一个乘法器输入端。然后将该乘法器输出的结果取反后与rx_yin(2)作加法运算，计算公式如下所示：

tx_yout(1)＝-rx_xin2^*×tx_xout+rx_yin(2)

选择器再将tx_xout*以及tx_xout输入到乘法器输入端，乘法器将输出的结果取反后与rx_yin(3)作加法运算，计算公式如下：

tx_yout(2)＝-tx_xout^*×tx_xout+rx_yin(3)

因此第二类方形PE包括两个乘法器、一个加法器和三个选择器。

第三类方形PE的结构框图如图7所示，该类方形PE共有4个输入，包括rx_xin1、rx_xin2，rx_xin3和rx_yin，第二类圆形PE的输出tx_xout_d与rx_xin1相连，第一类方形PE的输出tx_xout与rx_xin2相连，同时第二类方形处理单元的输出tx_xout与rx_xin3相连，并且rx_yin依次输入a₄₁,a₄₂,a₄₃,a₄₄。

其中输入端rx_xin1、rx_xin2与rx_xin3每M个时钟输入一个数据，rx_yin每M个时钟输入4个数据，分别为rx_yin(1)，rx_yin(2)，rx_yin(3)和rx_yin(4)。rx_xin1的第一个数据与rx_yin的第一个数据rx_yin(1)作乘法，然后将乘积作为tx_xout输出，计算公式如下：

tx_xout＝rx_yin(1)×rx_xin1

输出的同时，将tx_xout反馈回另一个到选择器的输入端，一个选择器的输入端输入全部tx_xout,另一个选择器的一个输入端输入tx_xout*，选择器将rx_xin2的共轭rx_xin2*以及tx_xout输入到另一个乘法器输入端。然后将该乘法器输出的结果取反后与输入rx_yin的第二个数据rx_yin(2)作加法运算，计算公式如下：

tx_yout(1)＝-rx_xin2^*×tx_xout+rx_yin(2)

依次类推，根据上述公式可以得到tx_yout(2)：

tx_yout(2)＝-rx_xin3^*×tx_xout+rx_yin(3)

在计算tx_yout(3)时，选择器再将tx_xout*以及tx_xout输入到乘法器输入端，乘法器输出的结果取反后与输入rx_yin的第三个数据rx_yin(3)作加法运算，计算公式如下：

tx_yout(3)＝-tx_xout^*×tx_xout+rx_yin(4)

所述第三类方形PE包括两个乘法器、一个加法器和三个选择器。

这样就完成了各个PE的结构设计。

步骤四：按照脉动阵列架构连接各个PE。

按照图1所示，将各个PE互联，形成满足要求的脉动阵列。在连接各个PE的时候，根据各个PE实际的输入和输出时序，添加必要的延迟单元，使得各个PE满足输入与输出的时序要求，实现数据的并行流水化处理。

步骤五：输出脉动阵列。

添加延迟单元后保证每一行PE的输出tx_xout作为脉动阵列在该行的有次序的输出，输出的顺序与PE所处的位置一致，从左到右按顺序依次进行，便于后续模块的接入或者数据观察。

按照上述步骤完成数据处理，便设计出完整的脉动阵列。

最后输出矩阵仍为(1，5，19)的定点形式，说明脉动阵列中的数据计算并没有出现错误。

在Verilog HDL的仿真中，选取仿真工具为Xilinx ISE 13.3及Modelsim SE6.6f，测试得到矩阵从完全输入到它的逆矩阵完全输出共需经历355个时钟。输入500组符合本发明实例的4阶正定Hermit随机复矩阵，将Verilog HDL仿真得到的Cholesky分解的结果与Matlab浮点仿真结果进行比较，计算每个矩阵中16个元素的平均相对误差，绘制出本发明方案下Verilog HDL定点仿真与Matlab浮点仿真的误差曲线，如图8所示。对这500个统计结果进一步取均值，得到平均误差为5.0545e-004。可以看到，采用本发明的设计方法，并且选择合适的定点方案，计算得到的结果具有比较高的精确度。

在Xilinx ISE中选取硬件平台XC7VX485T对本发明实例的工程进行综合，得到布局布线后的综合报告并将其占用的主要硬件资源列于表1中，其中，设计中的各乘法器以及除法器使用ISE中相应的IP核，且采用Mults方式实现。综合得到的最大时钟频率在280.513MHz。

表1本发明实例综合报告中主要硬件资源占有情况

资源类型	使用数量	可用数量	占有率
				Slice Registers	29,790	607,200	4％
Slice LUTs	14,572	303,600	4％
				bonded IOBs	385	600	64％
Block RAM/FIFO	39	1,030	3％
				DSP48E1s	72	2,800	2％

如上所述，便可以很好的实现本发明。

Claims (3)

1.正定Hermite矩阵Cholesky分解高速脉动阵列的设计方法，其特征在于，包括以下步骤：

(a)根据维数为M×M正定Hermite矩阵确定脉动阵列的维数以及构架；

(b)根据脉动阵列的架构对脉动阵列的处理单元进行分类；

(c)根据处理单元类别设计处理单元结构；

(d)根据脉动阵列的架构连接处理单元；

(e)输出脉动阵列；

所述步骤(a)中脉动阵列为M×M维的上三角单元处理阵列，且所述脉动阵列共有(1+M)M/2个处理单元；

所述步骤(b)中位于脉动阵列斜边上的处理单元被划分为圆形处理单元，其他单元被划分为方形处理单元；

所述步骤(c)中对圆形处理单元结构的设计具体包括以下步骤：

(c1)将脉动阵列上位于第M列、第M行的处理单元划分为第一类圆形处理单元，所述第一类圆形处理单元内置有一个开方器和一个寄存器，且包含一个输入端和一个输出端；

(c2)将脉动阵列斜边上的其他处理单元划分为第二类圆形处理单元，所述第二类圆形处理单元内置有一个开方器、一个除法器和两个寄存器，且包含一个输入端和两个输出端；

所述步骤(c)中对方形处理单元结构的设计包括以下步骤：

(I)将平行于圆形处理单元且与圆形处理单元相邻的处理单元划分为第一类方形处理单元，所述第一类方形处理单元由两个乘法器、一个加法器、一个选择器、两个移位寄存器和两个寄存器组成，且包含两个输入端和两个输出端；

(II)将平行于第一类方形处理单元且与第一类方形处理单元相邻的处理单元划分为第二类方形处理单元，所述第二类方形处理单元由两个乘法器、一个加法器、三个选择器两个移位寄存器和两个寄存器组成，且包含三个输入端和两个输出端；

(III)以此类推，将平行于第K-1类方形处理单元且与第K-1类方形处理单元相邻的处理单元划分为第K类方形处理单元，所述第K类方形处理单元由两个乘法器、一个加法器、三个选择器、两个移位寄存器和两个寄存器组成，其中第K类方形处理单元有K+1个输入端，两个输出端。

2.根据权利要求1所述的正定Hermite矩阵Cholesky分解高速脉动阵列的设计方法，其特征在于，所述K≤M，其中，K代表第K类方形处理单元，M代表维数为M×M正定Hermite矩阵中的维数。

3.根据权利要求2所述的正定Hermite矩阵Cholesky分解高速脉动阵列的设计方法，其特征在于，所述步骤(d)中在脉动阵列上的处理单元输出端还连接有延迟单元。