面向GPDSP的矩阵LU分解向量化计算的方法

摘要

一种面向GPDSP的矩阵LU分解向量化计算的方法，其步骤为：S1：根据GPDSP的体系结构特征确定最佳的LU分解的矩阵规模N值；S2：DSP核通过DMA从片外DDR存储器将要处理的矩阵数据传输到片内共享存储阵列中；S3：DSP核按照列选主元方法，计算第i列的主元值以及对应的列元素序号值；S4：根据列l

说明书

技术领域

本发明主要涉及通用计算数字信号处理器(General-Purpose Digital Signal Processor，简 称GPDSP)，特指一种适用于GPDSP的矩阵LU分解向量化计算的方法。

背景技术

稠密线性方程组求解是高性能计算和科学计算领域内最普遍的计算应用，而矩阵LU分 解(LU Factorization，LU)是求解稠密线性方程组最常用的一种方法，尤其是大规模稠密线 性方程组求解。高性能基准测试程序(High Performance Linpack，HPL)是TOP500最重要的 测试基准，HPL通过求解稠密线性代数方程组，以测试和评价高性能计算机系统的浮点性能。 HPL通过矩阵分块LU分解将大规模的矩阵分割为子块矩阵的LU分解计算，矩阵LU分解后 再通过两次三角方程组求解实现最终的方程组求解，其中矩阵LU分解的计算量占据整个HPL 计算量的95％以上。因此，优化矩阵LU分解的性能对提高HPL的效率具有非常重要的应用 价值。

在专利申请号为201310725118.6的文献(处于实审阶段)中提供了一种通用计算数字信 号处理器(General-Purpose Digital Signal Processor，简称GPDSP)，它包含CPU核单元和DSP 核单元，CPU核单元主要用于负责包括存储管理、文件控制、进程调度、中断管理任务在内 的通用事务管理以及提供对通用操作系统的完整支持；DSP核单元包含若干强大计算能力的 64位向量处理阵列，用于支持高密集运算任务的解算。传统的面向Cache结构的矩阵LU分 解方法不适合GPDSP的非Cache的向量阵列存储访存模式和向量处理阵列并发向量处理的体 系结构特征，难以发挥GPDSP向量计算优势。

发明内容

本发明要解决的技术问题就在于：针对现有技术存在的技术问题，本发明提供一种原理 简单、操作方便、能充分利用GPDSP的DSP核向量处理阵列的强大并行计算、高带宽向量 数据加载能力，显著提高DSP核计算访存比的面向GPDSP的矩阵LU分解向量化计算的方 法。

为解决上述技术问题，本发明采用以下技术方案：

一种面向GPDSP的矩阵LU分解向量化计算的方法，其步骤为：

S1：根据GPDSP的体系结构特征确定最佳的LU分解的矩阵规模N值；

S2：GPDSP的DSP核通过DMA从片外DDR存储器将需要处理的矩阵数据传输到片内 共享存储阵列中；

S3：GPDSP的DSP核按照列选主元方法，计算第i列的主元值以及对应的列元素序号值， 初始i＝0；

S4：GPDSP的DSP核根据上述计算得到的列l_i的列主元P_i以及对应的列元素序号值V_i对矩阵A内列主元P_i所在的行与列首元l_[i]所在的行进行交换；

S5：GPDSP的DSP核对列l_i按照公式l_i＝l_i/l_[i]进行列消元计算，更新列l_i；

S6：GPDSP的DSP核按照L’_i＝L’_i-l’_i*u_i进行更新矩阵panel的计算；

S7：判断i是否等于N-1，若不是，令i＝i+1，转步骤S3，若是转步骤S8；

S8：矩阵A的LU分解计算完毕；DSP核使用DMA将矩阵A从片内共享存储阵列传输 到片外DDR存储器的原存储位置。

作为本发明的进一步改进：在步骤S1中，设所处理的LU分解的矩阵规模为N×N阶非 奇异矩阵，每个矩阵元素的数据为u字节；上述矩阵规模的N值的确定方法是：设条件(1) 2*u*p*s*N≤q；(2)u*N*N≤r；按照满足条件(1)(2)取得的最大N值为所需要的N 值。

作为本发明的进一步改进：在步骤S2中，所述片内共享存储阵列中的矩阵标记为矩阵A， 其规模为N×N阶，首地址记为AA₀₀；沿左上至右下的对角线方向，将矩阵A划分为上三角 区域，标记为C，以及下三角区域，标记为D；将下三角矩阵D按照从左至右的方向按列分 别记为l₀、l₁、l₂…l_N-1，将每一列的在对角线上的数据，即列首元标记为l_[0]、l_[1]、l_[2]…l_[N-1]。

作为本发明的进一步改进：在所述步骤S5中，将片内共享存储阵列矩阵A中，列首元 l_[i]右下方的矩阵记为L’_i，矩阵A中第i行位于上三角C的部分记为u_i，列首元l_[i]右方包括 u_i部分的矩阵记为L_i，列l_i不包括列首元l_[i]的部分记为列l’_i。

作为本发明的进一步改进：所述步骤S3的详细流程为：

S3.1：DSP核通过DMA将列l_i从片内共享存储阵列中传输至片内向量阵列存储器；设初 始值为全0的向量寄存器Z₂和向量寄存器Z₅，以及初始值为{0,1,2,…,p-1}的向量寄存器Z₄；

S3.2：DSP核的向量处理阵列依次通过向量LOAD指令加载p个元素，存入向量寄存器 Z₀；

S3.3：DSP核的向量处理阵列对上述向量数据Z₀进行向量绝对值操作，结果存入向量寄 存器Z₁；

S3.4：DSP核的向量处理阵列对上述向量数据Z₁与Z₂的值进行向量比较操作，将比较 结果存入向量寄存器Z₃；

S3.5：根据上述比较结果值Z₃，DSP核的向量处理阵列使用向量MOV操作将Z₁中比Z₂数值更大的数据值更新向量寄存器Z₂；DSP核的向量处理阵列将数值更大的数据所对应存储 在Z₄中的列元素序号替换存入向量寄存器Z₅中；

S3.6：DSP核的向量处理阵列使用向量加法操作将向量寄存器Z₄的向量数值加p；

上述步骤中向量寄存器Z₄存储所处理列相对应的列元素序号；

重复步骤S3.1至步骤S3.6，直到完成该列的所有元素计算，最终得到p个列元素值以及 对应的p个列元素序号值；

S3.7：DSP核的向量处理阵列对上述p个列元素值进行数值大小的比较，计算出绝对值 最大的列元素值作为列l_i的列主元，标记为P_i，对应的列元素序号值标记为V_i。

作为本发明的进一步改进：所述步骤S4的详细流程为：

S4.1：根据列主元P_i的列元素序号值V_i，以及矩阵A在片内共享存储阵列中的首地址 AA₀₀，计算得出列主元P_i所在行的首地址AP_i＝AA₀₀+u*(V_i+i)*N；

S4.2：根据列首元l_[i]的列元素序号值，以及矩阵A在片内共享存储阵列中的首地址AA₀₀， 计算得出列首元l_[i]所在行的首地址Al_[i]＝AA₀₀+u*i*N；

S4.3：将片内向量阵列存储器划分为上下两片半区，半区的大小为q/2字节，DSP核使 用DMA将列主元P_i所在行共N个元素从片内共享存储阵列中传输至片内向量阵列存储器的 上半区；

S4.4：DSP核使用DMA将列首元l_[i]所在行共N个元素从片内共享存储阵列中传输至片 内向量阵列存储器的下半区；

S4.5：DSP核使用DMA将上述片内向量阵列存储器的下半区的N个元素数据传输至片 内共享存储阵列中以AP_i为首地址的位置；

S4.6：DSP核使用DMA将上述片内向量阵列存储器的上半区的N个元素数据传输至片 内共享存储阵列中以Al_[i]为首地址的位置。

作为本发明的进一步改进：所述步骤S5的详细流程为：

S5.1：DSP核使用DMA将列l_i除列首元素l_[i]以外的共N-i-1个元素数据从片内共享存储 阵列传输到片内向量阵列存储器；

S5.2：DSP核使用标量LOAD指令读取片内共享存储阵列内的列首元素l_[i]，存入标量 寄存器S₀；

S5.3：DSP核使用标量浮点除法指令计算1/l_[i]值，结果存入标量寄存器S₀，并使用标量 广播指令将S₀的数据广播至向量寄存器Z₀；

S5.4：DSP核使用向量LOAD指令依次加载上述已经传输至片内向量阵列存储器的该列 的其他元素数据，每次加载p个元素，存入向量寄存器Z₁；

S5.5：DSP核使用向量乘法指令将向量数据Z₁与向量数据Z₀相乘，结果存入向量寄存器 Z₁；

S5.6：DSP核使用向量STORE指令将Z₁中的向量数据结果存入片内向量阵列存储器中 的原位置；

重复步骤S5.4到步骤S5.6，直至列l_i的列消元计算完成；

S5.7：DSP核使用DMA将片内向量阵列存储器中更新后的列l_i传输至片内共享存储阵 列中的原位置。

作为本发明的进一步改进：将片内共享存储阵列矩阵A中，列首元l_[i]右下方的矩阵记为 L’_i，矩阵A中第i行位于上三角C的部分记为u_i，列首元l_[i]右方包括u_i部分的矩阵记为L_i， 列l_i不包括列首元l_[i]的部分记为列l’_i；所述步骤S6的详细流程为：

S6.1：DSP核使用DMA将列l’_i共N-i-1个元素从片内共享存储阵列传输到片内标量存储 器；

S6.2：将矩阵L_i按照规模v1×v2划分为w个分块矩阵，分块矩阵分别记为L_i,0，L_i,1，L_i,2，…， L_i，w-1，DSP核分w次循环，将每块分块矩阵传输至片内向量阵列存储器中进行计算，将计算 结果传回片内共享存储阵列的原存储位置；

上述矩阵L_i的分块规模v1，v2值的确定方法是：v1＝N-i，v2＝p；

上述循环次数w的确定方法是：w＝INT[(N-i-1)/p]，其中INT表示对方括号中的值向上取 整数；

将计算分为w次循环，第j+1次计算的分块矩阵记为L_i,j，矩阵L_i,j内包括矩阵L’_i的分块 矩阵记为L’_i,j，以及行u_i的部分记为u’_i，将DSP核的片内向量阵列存储器分为上下两片缓冲 区，缓冲区的大小为q/2字节，DSP核采用DMA双缓冲策略在片内向量阵列存储器和片内 共享存储阵列之间进行数据传输及计算。

作为本发明的进一步改进：所述步骤S6.2的具体流程为：

S6.2.1：DSP核使用DMA将分块矩阵L_i,j从片内共享存储阵列传输到片内向量存储器的 一个缓冲区中；

S6.2.2：DSP核启动DMA将下一块矩阵L_i,j+1从片内共享存储阵列传输到片内向量存储 器的另一个缓冲区中；

S6.2.3：DSP核在步骤S6.2.2启动DMA的同时，按照L’_i,j＝L’_i,j-l_i*u’_i计算更新矩阵L’_i,j；

上述步骤6.2.3按照如下步骤执行更新矩阵L’_i,j的计算；

S6.2.3.1：DSP核使用标量LOAD指获取片内标量存储器列l_i元素，存入标量寄存器S₀；

S6.2.3.2：DSP核使用标量广播指令将S₀数据广播至向量寄存器Z₀中；

S6.2.3.3：DSP核使用向量LOAD指令加载的p个u’_i向量数据，存入向量寄存器Z₁；

S6.2.3.4：DSP核使用向量LOAD指令加载的p个L’_i,j向量数据，存入向量寄存器Z₂；

S6.2.3.5：DSP核使用向量乘减指令执行Z₂-Z₁*Z₀；

S6.2.3.6：DSP核使用向量STORE指令将将计算结果Z₂存入片内向量存储器的原位置；

重复步骤S6.2.3.1至步骤S6.2.3.5，直至完成矩阵L’_i,j更新计算；

S6.2.4：DSP核启动DMA将更新后的矩阵L’_i,j从片内向量存储器的缓冲区中传输到片内 共享存储阵列的原位置；

重复步骤S6.2.1至步骤S6.2.4，直至完成矩阵L’_i更新计算。

与现有技术相比，本发明的优点在于：本发明的面向GPDSP的矩阵LU分解向量化计算 的方法，为针对GPDSP体系结构特征、高效的矩阵LU分解向量化计算方法，由DSP核的 向量处理阵列采用向量化方法计算串行的列选主元、行交换和列更新计算任务，由DSP核的 标量处理器和向量处理阵列紧密协同完成更新列panel的计算，充分发挥GPDSP的DSP核向 量处理阵列的强大并行计算和高带宽向量数据加载能力的优势，最终提高矩阵LU分解的计 算效率。

附图说明

图1是本发明在具体应用实例中所面向的GPDSP计算系统的简化存储模型示意图。

图2是本发明在具体应用实例中DSP核进行矩阵LU分解计算的流程示意图。

图3是本发明在具体应用实例中DSP核进行按列选取主元计算的流程示意图。

图4是本发明在具体应用实例中DSP核进行更新矩阵panel计算的流程示意图。

具体实施方式

以下将结合说明书附图和具体实施例对本发明做进一步详细说明。

如图1所示，为本发明在具体应用实例中所面向的GPDSP计算系统的简化存储模型示意 图。图中，GPDSP计算系统包含的CPU核单元、DSP核单元和片内共享存储阵列。DSP核 包含若干64位向量处理阵列，64位标量处理单元，片内标量存储器和片内向量阵列存储器。 片内共享存储阵列，用于CPU核与DSP核数据共享，并提供高带宽供数支持。片外DDR存 储器提供大容量的共享存储。

设GPDSP中DSP核的向量处理阵列计算单元数量为p个，每个计算单元包括s个MAC 部件(乘加部件)，DSP核的片内向量阵列存储器容量为q字节，GPDSP的片内共享存储阵 列容量为r字节，GPDSP的片外共享DDR存储容量为t字节。

如图2所示，本发明的面向GPDSP的矩阵LU分解向量化计算的方法，在应用实例中的 具体流程为：

S1：根据GPDSP的体系结构特征确定最佳的LU分解的矩阵规模N值。设所处理的LU 分解的矩阵规模为N×N阶非奇异矩阵，每个矩阵元素的数据为u字节。

上述矩阵规模的N值的确定方法是：设条件(1)2*u*p*s*N≤q；(2)u*N*N≤r； 按照满足条件(1)(2)取得的最大N值为所需要的N值。

S2：GPDSP的DSP核通过DMA从片外DDR存储器将需要处理的矩阵数据传输到片内 共享存储阵列中。

上述片内共享存储阵列中的矩阵标记为矩阵A，其规模为N×N阶，首地址记为AA₀₀。

沿左上至右下的对角线方向，将矩阵A划分为上三角区域(不包括对角线数据)，标记 为C，以及下三角区域(包括对角线数据)，标记为D。将下三角矩阵D按照从左至右的方向 按列分别记为l₀、l₁、l₂…l_N-1，将每一列的在对角线上的数据，即列首元标记为l_[0]、l_[1]、l_[2]…l_[N-1]。

S3：GPDSP的DSP核按照列选主元方法，计算第i列(初始i＝0)的主元值以及对应的 列元素序号值。

S4：GPDSP的DSP核根据上述计算得到的列l_i的列主元P_i以及对应的列元素序号值V_i对矩阵A内列主元P_i所在的行与列首元l_[i]所在的行进行交换。

S5：GPDSP的DSP核对列l_i按照公式l_i＝l_i/l_[i]进行列消元计算，更新列l_i。

将片内共享存储阵列矩阵A中，列首元l_[i]右下方的矩阵记为L’_i，矩阵A中第i行位于 上三角C的部分记为u_i，列首元l_[i]右方包括u_i部分的矩阵记为L_i，列l_i不包括列首元l_[i]的部 分记为列l’_i。

S6：GPDSP的DSP核按照L’_i＝L’_i-l’_i*u_i进行更新矩阵panel的计算。

S7：判断i是否等于N-1，若不是，令i＝i+1，转步骤S3，若是转步骤S8。

S8：矩阵A的LU分解计算完毕。DSP核使用DMA将矩阵A从片内共享存储阵列传输 到片外DDR存储器的原存储位置。

在具体应用实例中，如图3所示，上述步骤S3的详细流程为：

S3.1：DSP核通过DMA将列l_i从片内共享存储阵列中传输至片内向量阵列存储器。设初 始值为全0的向量寄存器Z₂和向量寄存器Z₅，以及初始值为{0,1,2,…,p-1}的向量寄存器Z₄。

S3.2：DSP核的向量处理阵列依次通过向量LOAD指令加载p个元素，存入向量寄存器 Z₀。

S3.3：DSP核的向量处理阵列对上述向量数据Z₀进行向量绝对值操作，结果存入向量寄 存器Z₁。

S3.4：DSP核的向量处理阵列对上述向量数据Z₁与Z₂的值进行向量比较操作，将比较 结果存入向量寄存器Z₃。

S3.5：根据上述比较结果值Z₃，DSP核的向量处理阵列使用向量MOV操作将Z₁中比Z₂ 数值更大的数据值更新向量寄存器Z₂；DSP核的向量处理阵列将数值更大的数据所对应存储 在Z₄中的列元素序号替换存入向量寄存器Z₅中。

S3.6：DSP核的向量处理阵列使用向量加法操作将向量寄存器Z₄的向量数值加p。

上述步骤中向量寄存器Z₄存储所处理列相对应的列元素序号。如列l_i的元素l_i,i、l_i,i+1、l_i+2,i对应的列元素序号为0,1,2，每次循环加p更新列元素序号。

重复步骤S3.1至步骤S3.6，直到完成该列的所有元素计算，最终得到p个列元素值以及 对应的p个列元素序号值。

S3.7：DSP核的向量处理阵列对上述p个列元素值进行数值大小的比较，计算出绝对值 最大的列元素值作为列l_i的列主元，标记为P_i，对应的列元素序号值标记为V_i。

在具体应用实例中，上述步骤S4的详细流程为：

S4.1：根据列主元P_i的列元素序号值V_i，以及矩阵A在片内共享存储阵列中的首地址 AA₀₀，计算得出列主元P_i所在行的首地址AP_i＝AA₀₀+u*(V_i+i)*N。

S4.2：根据列首元l_[i]的列元素序号值，以及矩阵A在片内共享存储阵列中的首地址AA₀₀， 计算得出列首元l_[i]所在行的首地址Al_[i]＝AA₀₀+u*i*N。

S4.3：将片内向量阵列存储器划分为上下两片半区，半区的大小为q/2字节，DSP核使 用DMA将列主元P_i所在行共N个元素从片内共享存储阵列中传输至片内向量阵列存储器的 上半区。

S4.4：DSP核使用DMA将列首元l_[i]所在行共N个元素从片内共享存储阵列中传输至片 内向量阵列存储器的下半区。

S4.5：DSP核使用DMA将上述片内向量阵列存储器的下半区的N个元素数据传输至片 内共享存储阵列中以AP_i为首地址的位置。

S4.6：DSP核使用DMA将上述片内向量阵列存储器的上半区的N个元素数据传输至片 内共享存储阵列中以Al_[i]为首地址的位置。

在具体应用实例中，上述步骤S5的详细流程为：

S5.1：DSP核使用DMA将列l_i除列首元素l_[i]以外的共N-i-1个元素数据从片内共享存储 阵列传输到片内向量阵列存储器。

S5.2：DSP核使用标量LOAD指令读取片内共享存储阵列内的列首元素l_[i]，存入标量 寄存器S₀。

S5.3：DSP核使用标量浮点除法指令计算1/l_[i]值，结果存入标量寄存器S₀，并使用标量 广播指令将S₀的数据广播至向量寄存器Z₀。

S5.4：DSP核使用向量LOAD指令依次加载上述已经传输至片内向量阵列存储器的该列 的其他元素数据，每次加载p个元素，存入向量寄存器Z₁。

S5.5：DSP核使用向量乘法指令将向量数据Z₁与向量数据Z₀相乘，结果存入向量寄存器 Z₁。

S5.6：DSP核使用向量STORE指令将Z₁中的向量数据结果存入片内向量阵列存储器中 的原位置。

重复步骤S5.4到步骤S5.6，直至列l_i的列消元计算完成。

S5.7：DSP核使用DMA将片内向量阵列存储器中更新后的列l_i传输至片内共享存储阵 列中的原位置。

在具体应用实例中，如图4所示，上述步骤S6的详细流程为：

将片内共享存储阵列矩阵A中，列首元l_[i]右下方的矩阵记为L’_i，矩阵A中第i行位于 上三角C的部分记为u_i，列首元l_[i]右方包括u_i部分的矩阵记为L_i，列l_i不包括列首元l_[i]的部 分记为列l’_i。

S6.1：DSP核使用DMA将列l’_i共N-i-1个元素从片内共享存储阵列传输到片内标量存储 器。

S6.2：将矩阵L_i按照规模v1×v2划分为w个分块矩阵，分块矩阵分别记为L_i,0，L_i,1，L_i,2，…， L_i，_w-1，DSP核分w次循环，将每块分块矩阵传输至片内向量阵列存储器中进行计算，将计算 结果传回片内共享存储阵列的原存储位置。

上述矩阵L_i的分块规模v1，v2值的确定方法是：v1＝N-i，v2＝p。

上述循环次数w的确定方法是：w＝INT[(N-i-1)/p]，其中INT表示对方括号中的值向上取 整数。

将计算分为w次循环，第j+1次计算的分块矩阵记为L_i,j，矩阵L_i,j内包括矩阵L’_i的分块 矩阵记为L’_i,j，以及行u_i的部分记为u’_i，将DSP核的片内向量阵列存储器分为上下两片缓冲 区，缓冲区的大小为q/2字节，DSP核采用DMA双缓冲策略在片内向量阵列存储器和片内 共享存储阵列之间进行数据传输及计算。

S6.2.1：DSP核使用DMA将分块矩阵L_i,j从片内共享存储阵列传输到片内向量存储器的 一个缓冲区中。

S6.2.2：DSP核启动DMA将下一块矩阵L_i,j+1从片内共享存储阵列传输到片内向量存储 器的另一个缓冲区中。

S6.2.3：DSP核在步骤S6.2.2启动DMA的同时，按照L’_i,j＝L’_i,j-l_i*u’_i计算更新矩阵L’_i,j。

上述步骤6.2.3按照如下步骤执行更新矩阵L’_i,j的计算。

S6.2.3.1：DSP核使用标量LOAD指获取片内标量存储器列l_i元素，存入标量寄存器S₀。

S6.2.3.2：DSP核使用标量广播指令将S₀数据广播至向量寄存器Z₀中。

S6.2.3.3：DSP核使用向量LOAD指令加载的p个u’_i向量数据，存入向量寄存器Z₁。

S6.2.3.4：DSP核使用向量LOAD指令加载的p个L’_i,j向量数据，存入向量寄存器Z₂。

S6.2.3.5：DSP核使用向量乘减指令执行Z₂-Z₁*Z₀。

S6.2.3.6：DSP核使用向量STORE指令将将计算结果Z₂存入片内向量存储器的原位置。

重复步骤S6.2.3.1至步骤S6.2.3.5，直至完成矩阵L’_i,j更新计算。

S6.2.4：DSP核启动DMA将更新后的矩阵L’_i,j从片内向量存储器的缓冲区中传输到片内 共享存储阵列的原位置。

重复步骤S6.2.1至步骤S6.2.4，直至完成矩阵L’_i更新计算。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于 本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术 人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。