, 吴艳霞 , 张博为 哈尔滨工程大学 计算机科学与技术学院, 哈尔滨 150001
收稿日期: 2016-04-11
基金项目: 国家自然科学基金面上项目（61003036）；计算机体系结构国家重点实验室开放课题（CARCH201301）；中央高校基本科研业务经费专项基金（HEUCF100606）
作者简介: 刘书勇(1978-),男,讲师
通讯作者: 林俊宇,讲师,E-mail:linjunyu@hrbeu.edu.cn
摘要：矩阵运算是高性能计算中核心问题之一，矩阵分解是提高矩阵运算并行性的重要途径，飞速发展的FPGA为并行运算结构提供了有力的环境支持。该文基于子矩阵更新同一化算法实现了Cholesky分解，基于FPGA设计了相应的并行结构。实验结果表明：与通用处理器的软件实现相比，本文实现的Cholesky分解的FPGA并行结果在核心计算性能上可以取得10倍以上的加速比，该算法针对矩阵三角化计算过程具有更高的数据和流水并行性。
关键词： 矩阵三角化分解 Cholesky分解 并行结构 现场可编程门阵列
Cholesky decomposition and parallel structure design based on matrix triangularization decomposition
LIU Shuyong, LIN Junyu
, WU Yanxia, ZHANG BoweiCollage of Computer Science and Technology, Harbin Engineering University, Harbin 150001, China
Abstract:Matrix computing is one of the core problems in high performance computing with matrix decomposition being an important way to improve the parallelism of matrix computations. FPGA gives a powerful environment for parallel computing. This study uses Cholesky decomposition based on a hardware-adaptive parallel sub-matrix identity updating algorithm. Its parallel structure is based on FPGA. Tests show that this structure achieves more than 10 fold speedup compared to general-purpose processors in terms of the kernel computational speed because the algorithm has better data-parallelism and pipeline-parallelism during matrix triangularization.
Key words: matrix triangularization decompositionCholesky decompositionparallel structurefield programmable gate array
当前，对高性能矩阵三角化分解的研究主要从通用计算[1]、脉动阵列[2]和FPGA实现[3]3个方向开展。在高性能通用计算领域内，从计算特征及应用领域考虑，对矩阵三角化分解的研究主要在基于单指令流多数据流(SIMD)或多指令流多数据流(MIMD)的向量机、共享存储的多处理机等技术[4-10]的并行计算机体系结构上开展。一些研究项目在单处理器以及通用计算图形处理器(GPGPU)上也有良好的表现。该领域的研究通过向程序设计者提供具有标准接口和较好移植性的数值计算模块，以在目标体系结构上实现算法的高性能计算。许多高层编程语义数值计算库如基于C++的cvmlib、英特尔数学核心函数库(Intel MKL)和矩阵计算模板库(matrix template library,MTL)等则通过调用1个或若干基础数值计算子程序库来实现。在业界已有的研究成果中，以BLAS、LAPACK和LINPACK等[11-12]为代表的一系列数值线性代数软件包最具有代表性。使用线性或二维脉动阵列机实现三角化分解的研究集中在20世纪80年代和90年代初，往往针对某些特定科学工程领域的应用。该结构是VLSI设计中一种经典体系结构，通过消除阵列中的广播路径实现三角化分解中流水并行的最大化，往往占用较大的片上面积，产生较大的功耗。近年来，更多的研究者在固定面积的FPGA上设计具有较高性能的矩阵三角化分解并行结构。归纳起来，这些设计分别选择从单个处理单元(processing element,PE)设计、存储层次设计、阵列结构设计等方面开展优化，实现FPGA矩阵三角化分解的研究以及针对软件数值线性代数软件包的FPGA设计[13-15]。从实现结构角度来说，它们大多是针对某种具体矩阵分解的一维线性阵列及二维多边形阵列以不同并行算法加以实现[16]。这些设计在科学及工程应用领域已经得到广泛应用。但是，对FPGA矩阵三角化分解的研究仍有可提升的空间。从矩阵计算类属的角度考虑，专用的矩阵分解并行算法及阵列结构很难直接复用到其他矩阵分解的实现中，而以通用化为设计目标的结构则在计算性能上无法与专用阵列结构相比。此外，已有的一些矩阵分解细粒度并行算法往往需要对算法本身进行较大的改动，通过增加硬件操作原语来实现算法指令级的并行[17-18]。在设计一类具有共同计算特征的矩阵分解应用中，对计算过程和算法本身细粒度并行的分析和优化成为提高计算阵列性能的关键。矩阵三角化分解过程中均含有三角化计算，这一数值计算过程同样具有规律性较强的数据排列与流动特征，非常适合进行较细粒度的FPGA并行优化设计。
本文针对矩阵三角化分解中共有的三角化计算迭代过程，通过分析其中子矩阵更新过程的线性规律，提出一种具有一般性的子矩阵更新同一化并行算法，进而提出基于该算法的计算复杂度较低的矩阵三角化分解并行结构，针对Cholesky矩阵三角化分解在并行结构上的高性能实现及优化方法开展研究。
1 研究背景 本文提到的矩阵和向量均在实数范围内，矩阵皆为方阵。
单个矩阵的三角化分解将矩阵A分解为正交矩阵与三角矩阵之积，或分解为一个上三角矩阵与一个下三角矩阵之积。三角化过程可以总结为： 给定n×p(p≤n)维的矩阵X，根据式(1)求n×n维的非奇异矩阵M以及p×p维的上三角矩阵R。
| $MX=\left[ \begin{align} & R \\ & 0 \\ \end{align} \right].$ | (1) |
2 Cholesky分解及FPGA并行结构设计 Cholesky分解将n×n维对称正定矩阵A(其中元素为aij,1≤i,j≤n)分解成A=LLT的形式，L(其中元素为lij,1≤j≤i≤n)为对角线元素为正数的下三角矩阵。L与LT如式(2)所示。
| $\begin{align} & L=\left[ \begin{matrix} {{l}_{11}} & 0 & \cdots & 0 & \cdots & 0 \\ {{l}_{21}} & {{l}_{22}} & \cdots & 0 & \cdots & 0 \\ \vdots & \vdots & {} & \vdots & {} & \vdots \\ {{l}_{k1}} & {{l}_{k2}} & \cdots & {{l}_{kk}} & \cdots & {} \\ \vdots & \vdots & {} & \vdots & {} & \vdots \\ {{l}_{n1}} & {{l}_{n2}} & \cdots & {{l}_{nk}} & \cdots & {{l}_{nn}} \\\end{matrix} \right], \\ & {{L}^{T}}=\left[ \begin{matrix} {{l}_{11}} & {{l}_{21}} & \cdots & {{l}_{k1}} & \cdots & {{l}_{n1}} \\ 0 & {{l}_{22}} & \cdots & {{l}_{k2}} & \cdots & {{l}_{n2}} \\ \vdots & \vdots & {} & \vdots & {} & \vdots \\ 0 & 0 & \cdots & {{l}_{kk}} & \cdots & {{l}_{nk}} \\ \vdots & \vdots & {} & \vdots & {} & \vdots \\ 0 & 0 & \cdots & 0 & \cdots & {{l}_{nn}} \\\end{matrix} \right]. \\ \end{align}$ | (2) |
| ${{l}_{ij}}=\left\{ \begin{matrix} {{a}_{ij}}/{{l}_{jj}}, & italic>j=1; \\ \left( {{a}_{ij}}-\sum\limits_{k=1}^{j-1}{{{l}_{ik}}{{l}_{jk}}} \right)/{{l}_{jj}}, & italic>j>1; \\\end{matrix} \right.~$ | (3) |
| ${{l}_{ii}}=\sqrt{{{a}_{ii}}-\sum\limits_{k=1}^{j-1}{l_{ik}^{2}}}.$ | (4) |
 |
| 图 1 LT-SC的计算过程 |
| 图选项 |
因此，LT-SC基于列k的1次子矩阵更新亦可以表示为式(5)和(6)。
| $\begin{align} & {{L}^{T}}_{(1,1)}=\left[ \begin{matrix} {{u}_{11}} & {{u}_{12}} & {{u}_{13}} & {{u}_{14}} \\ \times & a_{22}^{'} & a_{23}^{'} & a_{24}^{'} \\ \times & \times & a_{33}^{'} & a_{34}^{'} \\ \times & \times & \times & a_{44}^{'} \\\end{matrix} \right] \\ & ,{{\left[ \begin{matrix} {{u}_{11}} \\ {{u}_{12}} \\ {{u}_{13}} \\ {{u}_{14}} \\\end{matrix} \right]}^{T}}=\left[ \begin{matrix} \sqrt{{{a}_{11}}} \\ {{a}_{12}}/\sqrt{{{a}_{11}}} \\ {{a}_{13}}/\sqrt{{{a}_{11}}} \\ {{a}_{14}}/\sqrt{{{a}_{11}}} \\\end{matrix} \right], \\ & {{A}^{'}}=\left[ \begin{matrix} {{a}_{11}} & {{a}_{12}} & {{a}_{13}} & {{a}_{14}} \\ \times & a_{22}^{'} & a_{23}^{'} & a_{24}^{'} \\ \times & \times & a_{33}^{'} & a_{34}^{'} \\ \times & \times & \times & a_{44}^{'} \\\end{matrix} \right]= \\ & {{T}_{1,1}}\left[ \begin{matrix} {{a}_{11}} & {{a}_{12}} & {{a}_{13}} & {{a}_{14}} \\ {{a}_{21}} & {{a}_{22}} & {{a}_{23}} & {{a}_{24}} \\ {{a}_{31}} & \times & {{a}_{33}} & {{a}_{34}} \\ {{a}_{41}} & \times & \times & {{a}_{44}} \\\end{matrix} \right] \\ \end{align}$ | (5) |
| ${T_{1,1}} = \left[ {\matrix{ 1 & 0 & 0 & 0 \cr { - {a_{21}}/{a_{11}}} & 1 & 0 & 0 \cr { - {a_{31}}/{a_{11}}} & 0 & 1 & 0 \cr { - {a_{41}}/{a_{11}}} & 0 & 0 & 1 \cr } } \right].$ | (6) |
| $P=\left[ \begin{matrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \\\end{matrix} \right]$ | (7) |
| $\left\{ \begin{align} & {{A}^{(t)}}=P({{T}_{1,1}}{{A}^{(t-1)}})P,0 ＜t\le n. \\ & {{A}^{(0)}}=A \\ \end{align} \right.$ | (8) |
| $z_{~i-1,j-1}^{(t)}={{g}^{(t)}}\left( i,j,k \right)=z_{i,j}^{(t-1)}-\left( \text{ }z_{i,1}^{(t-1)}/\text{ }z_{1,1}^{(t-1)} \right)z_{1,j}^{(t-1)},$ |
则使用子矩阵同一化并行算法后，子矩阵的更新过程可以表示为式(9)。
| ${T_d}{A_{4 \times 4}} = {T_d}\left[ {\matrix{ {{a_{11}}} & {{a_{12}}} & {{a_{13}}} & {{a_{14}}} \cr \times & {{a_{22}}} & {{a_{23}}} & {{a_{24}}} \cr \times & \times & {{a_{33}}} & {{a_{34}}} \cr \times & \times & \times & {{a_{44}}} \cr } } \right],$ | (9a) |
| $\eqalign{ & {T_d} = \left[ {\matrix{ { - {a_{12}}/{a_{11}}} & 1 & 0 & 0 \cr { - {a_{13}}/{a_{11}}} & 0 & 1 & 1 \cr { - {a_{14}}/{a_{11}}} & 0 & 1 & 1 \cr 1 & 0 & 1 & 1 \cr } } \right] = \cr & p\left[ {\matrix{ 1 & 0 & 0 & 0 \cr { - {a_{12}}/{a_{11}}} & 1 & 0 & 0 \cr { - {a_{13}}/{a_{11}}} & 0 & 1 & 0 \cr { - {a_{14}}/{a_{11}}} & 0 & 0 & 1 \cr } } \right]. \cr} $ | (9b) |
 |
| 图 2 LT-SC |
| 图选项 |
 |
| 图 3 LT-SC的阵列结构与计算时空图 |
| 图选项 |
 |
| 图 4 LT-SC并行结构PE的配置 |
| 图选项 |
 |
| 图 5 LT-SC总体并行结构 |
| 图选项 |
3 实验结果本文提出的设计采用VHDL语言实现。使用Mentor Graphic Modelsim v6.5仿真验证工具进行仿真后在Xilinx ISE 11.1工具中完成综合，并给出资源占用结果； 对设计的性能进行估计与分析； 最后将设计加载到FPGA上进行实际运行与功能验证。所使用的FPGA分别为xc5vlx330t和xc5vlx110t。
3.1 FPGA资源占用本文涉及的浮点计算如加法、减法、乘法、除法以及开平方等均使用开源算术浮点库libhdlfltp的单精度流水浮点计算部件实现，遵照IEEE754标准的浮点数格式。
表 1给出LT-SC并行结构中PE单元的综合结果，包括BcPE-1、IcPE、BcPE-2和rPE等4种PE单元的资源使用和最大频率。可以看到，因内部没有配置浮点运算部件，rPE的资源占用非常少。
表 1 LT-SC并行结构中PE单元的综合结果
| 单元 | 寄存器个数 | 查找表个数 | 查找表触发器个数 | DSP48E个数 | 最大频率/MHz |
| BcPE-1 | 1 156 | 1 212 | 1 228 | 0 | 188.56 |
| IcPE | 338 | 847 | 868 | 2 | 169.21 |
| BcPE-2 | 2 034 | 1 942 | 2 054 | 0 | 150.76 |
| rPE | 34 | 40 | 42 | 0 | 550.01 |
表选项
3.2 性能比较用于比较的台式计算机配置为主频2.33 GHz 英特尔酷睿2处理器、4 GB主存及CentOS 5.5操作系统。C语言编译器为gcc 4.3，优化选项设为O2。图 6给出LT-SC阵列并行结构与通用处理器上软件实现相比所获得的加速比和最大频率随矩阵规模(阵列规模)增加的变化。可以看到，在各个矩阵规模上的计算中LT-SC阵列并行结构均可以取得10倍以上的加速比，即使最大频率在降低，加速比也随着矩阵规模的增加逐步提升。需要指出的是，将矩阵规模从10×10维增加到18×18维过程中，随着IcPE数量的大量增加，FPGA上的DSP48E Slice资源全部消耗，故最大频率迅速下降。而矩阵规模接近32×32维时，随着FPGA上LUT Flip-Flop pairs资源的耗尽，最大频率大幅度下降。
 |
| 图 6 LT-SC阵列并行结构与通用处理器软件实现的加速比和最大频率随矩阵规模增加的变化 |
| 图选项 |
4 结 论本文基于子矩阵更新同一化算法，针对Cholesky分解的求解特征，提出以求解矩阵L的转置阵LT的LT-SC算法，并在分析LT-SC线性代数特征的基础上，给出了LT-SC的并行结构实现方案。在并行结构的具体实现中，采用FPGA协处理器浮点数学库libhdlfltp支持阵列PE的浮点流水运算。理论分析表明，该算法针对矩阵三角化计算过程具有更高的数据并行性与流水并行性； 实验结果表明，针对Cholesky分解，基于LT-SC算法在FPGA上的实现在核心计算性能上优于在通用处理器上的软件实现，可达到10倍以上的加速比。但作为实际运行硬件的可重构计算验证平台中的通信开销所占据整体运行时间比例过大，已经成为整体计算加速的瓶颈，下一步应加以改进。
参考文献
| [1] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Satchidanand G H, Sotirios G Z. FPGA implementation of a Cholesky algorithm for a shared-memory multiprocessor architecture[J]. Parallel Algorithms and Applications, 2004, 19(4) : 211–226.DOI:10.1080/10637190412331279957 |
| [2] | Journal of Central South University(Science and Technology), 41(2):649-654.--> Kung H T, Leiserson C E. Systolic arrays technical report [R]. Pitsburg, PA, USA: Carneige Mellon University, 1978. |
| [3] | Journal of Central South University(Science and Technology), 41(2):649-654.-->YANG Depeng, Gregory D P, LI Husheng. Compressed sensing and Cholesky decomposition on FPGAs and GPUs[J]. Parallel Computing, 2012, 38(8) : 421–437.DOI:10.1016/j.parco.2012.03.001 |
| [4] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Alfredo B, Julien L, Jakub K, et al. A class of parallel tiled linear algebra algorithms for multicore architectures[J]. Parallel Computing, 2009, 35(1) : 38–53.DOI:10.1016/j.parco.2008.10.002 |
| [5] | Journal of Central South University(Science and Technology), 41(2):649-654.--> DOU Yong, Vassiliadis S, Kuzmanov G, et al. 64-bit floating-point FPGA matrix multiplication [C]//Proceedings of the 13th ACM/SIGDA International Symposium on Field Programmable Gate Arrays(FPGA’05), New York, NY, USA: ACM, 2005: 86-95. http://cn.bing.com/academic/profile?id=2061624656&encoded=0&v=paper_preview&mkt=zh-cn |
| [6] | Journal of Central South University(Science and Technology), 41(2):649-654.--> Vasily V, James W D. LU, QR and Cholesky factorizations using vector capabilities of GPUs [R]. Berkeley, CA, USA: EECS Department, University of California, Berkeley, 2008. http://cn.bing.com/academic/profile?id=2187691968&encoded=0&v=paper_preview&mkt=zh-cn |
| [7] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Oleg M, Volodymyr L, Anatoli S, et al. Parallel implementation of Cholesky LLT-algorithm in FPGA-based processor[J]. Parallel Processing and Applied Mathematics, 2008, 49(67) : 137–147. |
| [8] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Eunice E S, CHU Peiyun. Efficient and optimal parallel algorithms for Cholesky decomposition[J]. Journal of Mathematical Modeling and Algorithms, 2003, 2(3) : 217–234.DOI:10.1023/B:JMMA.0000015832.41014.ed |
| [9] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Aatonio R, George A. A high throughput FPGA-based floating point conjugate gradient implementation for dense matrices[J]. ACM Transactions on Reconfigurable Technology and Systems, 2010, 3(1) : 1–19. |
| [10] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Badia J M, Vidal A M. Parallel algorithms to computer the eigenvalues and eigenvectors of symmetric Toeplitz matrices[J]. Parallel Algorithms and Applications, 1998, 13(1) : 75–93.DOI:10.1080/01495739808947361 |
| [11] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Anderson E, BAI Z, Bischof C, et al. LAPACK Users Guide[M].Philadelphia, PA, USA: The society of industrial and applied mathematics, 1999. |
| [12] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Blackford L S, Choi J, Cleary A, et al. ScaLAPACK Users Guide[M].Philadelphia, PA, USA: The Society of Industrial and Applied Mathematics, 1997. |
| [13] | Journal of Central South University(Science and Technology), 41(2):649-654.-->WU Guiming, DOU Yong, SUN Junqing, et al. A high performance and memory efficient LU decomposer on FPGAs[J]. Computers, IEEE Transactions, 2012, 61(3) : 366–378.DOI:10.1109/TC.2010.278 |
| [14] | Journal of Central South University(Science and Technology), 41(2):649-654.--> WU Guiming, DOU Yong, Gregory D P. Blocking LU decomposition for FPGAs [C]//Proceedings of the 18th IEEE Symposium on Field-Programmable Custom Computing Machines(FCCM’10). Charlotte, NC,USA: IEEE, 2010, 109-112. |
| [15] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Haridas S G, Sotirios G Z. FPGA implementation of a Cholesky algorithm for a shared memory multiprocessor architecture[J]. International Journal of Parallel, Emergent and Distributed Systems, 2004, 19(4) : 211–226. |
| [16] | Journal of Central South University(Science and Technology), 41(2):649-654.-->郭磊, 唐玉华, 周杰, 等. 基于FPGA的Cholesky分解细粒度并行结构与实现[J]. 计算机研究与发展, 48 (Suppl), 2011, 48(supp) : 258–265.GUO Lei, TANG Yuhua, ZHOU Jie, et al. Fine-grained architecture and implementation for Cholesky decomposition on FPGA[J]. Journal of Computer Research and Development, 2011, 48(supp) : 258–265.(in Chinese) |
| [17] | Journal of Central South University(Science and Technology), 41(2):649-654.-->WANG Xiaojun, Leeser M. A truly two-dimensional systolic array FPGA implementation of QR decomposition[J]. ACM Transactions on Embedded Computing Systems (TECS), 2009, 9(1) : 17–18. |
| [18] | Journal of Central South University(Science and Technology), 41(2):649-654.-->邬贵明, 窦勇, 王淼. Cholesky分解细粒度并行算法[J]. 计算机工程与科学, 2010, 32(9) : 102–106.WU Guiming, DOU Yong, WANG Miao. A fine-grained parallel algorithm for the Cholesky decomposition[J]. Computer Engineering & Science, 2010, 32(9) : 102–106.(in Chinese) |
| [19] | Journal of Central South University(Science and Technology), 41(2):649-654.-->刘书勇, 吴艳霞, 张博文, 等. 基于可重构计算系统的矩阵三角化分解硬件并行结构研究[J]. 电子学报, 43(8), 2015, 43(8) : 1642–1650.LIU Shuyong, WU Yanxia, ZHANG Bowen, et al. Research of parallel hardware architecture for matrix triangularization decomposition based on reconfigurable computing system[J]. ACTA Electronica Sinica, 2015, 43(8) : 1642–1650.(in Chinese) |
