基于GPU的稀疏矩阵LU分解方法.pdf

《基于GPU的稀疏矩阵LU分解方法.pdf》由会员分享，可在线阅读，更多相关《基于GPU的稀疏矩阵LU分解方法.pdf（18页珍藏版）》请在专利查询网上搜索。

1、(10)申请公布号 CN 103399841 A (43)申请公布日 2013.11.20 CN 103399841 A *CN103399841A* (21)申请号 201310329479.9 (22)申请日 2013.07.31 G06F 17/16(2006.01) (71)申请人 清华大学 地址 100084 北京市海淀区 100084-82 信箱 (72)发明人 任令 陈晓明 汪玉 杨华中 (74)专利代理机构 北京清亦华知识产权代理事 务所 ( 普通合伙 ) 11201 代理人 张大威 (54) 发明名称 基于 GPU 的稀疏矩阵 LU 分解方法 (57) 摘要 本发明提出一种基。

2、于GPU的稀疏矩阵LU分解 方法， 包括以下步骤 ： A ： 在CPU上建立待仿真的电 路网单的电路矩阵， 并对该电路矩阵进行预处理 和稀疏化 ； B ： 根据预处理结果选择对稀疏矩阵 LU 进行处理的计算平台 ； C ： 如果使用的计算平台为 GPU， 则通过 GPU 对稀疏矩阵 LU 进行分解 ； D ： 通过 CPU 对稀疏矩阵 LU 中非零元的值进行回代求解， 以完成对稀疏矩阵LU的分解。 本发明的实施例可 根据分解中浮点运算次数的多少自动选择计算平 台， 能够合理分配 CPU 和 GPU 上的任务， 保证 GPU 上众多线程间的协同运行， 从而提高运算速度， 减 少电路仿真时间。 (。

3、51)Int.Cl. 权利要求书 2 页 说明书 10 页 附图 5 页 (19)中华人民共和国国家知识产权局 (12)发明专利申请 权利要求书2页 说明书10页 附图5页 (10)申请公布号 CN 103399841 A CN 103399841 A *CN103399841A* 1/2 页 2 1. 一种基于 GPU 的稀疏矩阵 LU 分解方法， 其特征在于， 包括以下步骤 ： A ： 在 CPU 上建立待仿真的电路网单的电路矩阵， 并对所述电路矩阵进行预处理和稀疏 化 ； B ： 根据预处理结果选择对所述稀疏矩阵 LU 进行处理的计算平台 ； C ： 如果对所述稀疏矩阵 LU 进行处理所。

4、使用的计算平台为 GPU， 则通过所述 GPU 对所述 稀疏矩阵 LU 进行分解， 包括 ： C1 ： 将所述电路矩阵的非零元位置、 值以及所述稀疏矩阵的 LU 因子的非零元位置输入 所述 GPU ； C2 ： 确定群并行模式和流水线并行模式的分界点 ； C3 ： 依次计算属于所述群并行模式的各层直至所有层均计算完成 ； C4 ： 通过所述流水线并行方式计算剩余各层 ； C5 ： 将所述稀疏矩阵 LU 中非零元的值写回到所述 CPU ； D ： 通过所述 CPU 对所述稀疏矩阵 LU 中非零元的值进行回代求解， 以完成对稀疏矩阵 LU 的分解。 2. 如权利要求 1 所述的基于 GPU 的稀疏。

5、矩阵 LU 分解方法， 其特征在于， 所述对所述电 路矩阵进行预处理， 进一步包括 ： 对数值进行处理以提高所述数值的稳定性 ； 根据近似最小自由度算法减少分解过程中的非零元写入 ； 根据部分选主元的数值分解计算所述矩阵 LU 的非零元结构。 3. 如权利要求 2 所述的基于 GPU 的稀疏矩阵 LU 分解方法， 其特征在于， 根据 HSL_MC64 方法对所述数值进行处理以提高所述数值的稳定性。 4. 如权利要求 1 所述的基于 GPU 的稀疏矩阵 LU 分解方法， 其特征在于， 所述根据预处 理结果选择对所述稀疏矩阵 LU 进行处理的计算平台， 进一步包括 ： 对所述矩阵 LU 中每行的非。

6、零元按列号进行排序 ； 建立消去图， 并对所述消去图进行分层 ； 根据浮点运算次数确定要使用的平台 ； 如果所述浮点运算次数大于或等于预设次数， 则使用 GPU 平台 ； 如果所述浮点运算次数小于所述预设次数， 则使用 CPU 平台。 5. 如权利要求 4 所述的基于 GPU 的稀疏矩阵 LU 分解方法， 其特征在于， 所述建立消去 图， 并对所述消去图进行分层， 进一步包括 ： 建立消去图， 其中， 所述消去图中包括节点， 且各节点表示各行。 对所述消去图进行分层， 其中， 所述节点的层号为所述节点所依赖的行所处的最大层 号加 1。 6. 如权利要求 5 所述的基于 GPU 的稀疏矩阵 LU。

7、 分解方法， 其特征在于， 所述对所述消 去图进行分层之后， 还包括 ： 依次处理所述消去图的每一层， 所述每一层内部的各行可按照任意顺序进行计算。其 中， 处理所述消去图的方式包括 ： 群并行模式和流水线并行模式， 所述群并行模式指所述 每一层有足够多的行被并行分解 ； 所述流水线并行模式指所述每一层只有很少的行， 且来 权 利 要 求 书 CN 103399841 A 2 2/2 页 3 自不同层的行重叠执行分解。 7. 如权利要求 1 所述的基于 GPU 的稀疏矩阵 LU 分解方法， 其特征在于， 所述通过所述 流水线并行方式计算剩余各层， 进一步包括 ： 如果检测到未计算完成的行， 则。

8、设为第 k 行， 并计算所述第 k 行 ； 依次检测所述第 k 行所依赖的行， 并设为第 j 行， 并判断所述第 j 行是否计算完成 ； 如果所述第 j 行计算完成， 则从第 K 行中减去第 j 行乘以其相应系数 ； 如果第 j 行未完成， 则等待预设时间后重新检测第 j 行是否计算完成 ； 依次确认所述第 k 行所依赖的每一行均计算完成后， 从所述第 k 行中减去该行乘以其 相应系数 ； 将第 k 行的计算结果写回到所述 GPU 中， 并判定第 k 行已计算完成， 并进一步处理下一 个未完成的行， 直至所有行均已完成计算。 权 利 要 求 书 CN 103399841 A 3 1/10 页 。

9、4 基于 GPU 的稀疏矩阵 LU 分解方法 技术领域 0001 本发明涉及电子设计自动化与并行计算领域， 特别涉及一种基于 GPU 的稀疏矩阵 LU 分解方法。 背景技术 0002 稀疏矩阵 LU 分解是线性代数中的一项基本操作， 在电路仿真、 结构力学、 经济建 模等许多领域有广泛应用。在电路仿真中， 电路被表示为矩阵， 电路矩阵极其稀疏， 通常每 行只有约 5 个非零元， 这是由在电路中 （除少量电源、 地、 时钟节点外） 一个节点不能连接太 多元件这一原因决定的。求解电路的过程， 也就是利用稀疏矩阵 LU 分解求解稀疏矩阵方程 Ax=b的过程。 在一次电路仿真中， 稀疏矩阵LU分解位于。

10、牛顿拉夫森迭代和瞬态迭代两层 循环之中， 需要迭代地进行很多次， 这一步骤是当前最常用的电路仿真 SPICE(Simulation Program with Integrated Circuit Emphasis) 的主要瓶颈。目前对超大规模集成电路的 仿真的时间可能长达数周甚至数月。因此稀疏矩阵 LU 分解的并行化是一个急需解决的问 题。 0003 矩阵 LU 分解的目标是找到上三角矩阵 L 和下三角矩阵 U 使得 0004 0005 上述方程解为 0006 0007 0008 上式可以写成 0009 说 明 书 CN 103399841 A 4 2/10 页 5 0010 lij=xij,。

11、j i 0011 0012 当矩阵 A 是稀疏矩阵时， 其 LU 因子矩阵 L 和 U 往往也是稀疏矩阵， 此时， 当且仅当 lik 0 时， 需要从 xi 中减去 uk。 0013 LU 分解完成后， 求解 Ax=b 的问题即变成求解 Ly=b 和 Ux=y 两个三角方程的问题， 这一步称为回代求解。 0014 理论上， 将稀疏矩阵当作普通稠密矩阵进行LU分解， 也可得到正确结果的。 加州大 学伯克利分校 （University of California, Berkeley） 、 田纳西大学诺克斯维尔分校 （The University of Tennessee, Knoxville） 的。

12、研究人员已经能在 GPU 上高效的实现稠密矩阵 的LU分解， 在GTX 280上达到388十亿次每秒(Gflop/s,Giga Floating point OPerations per second)，相 关 研 究 可 参 见 “V.Volkov and J.Demmel,“LU,QR and Cholesky factorizations using vector capabilities of gpus,” EECS Department,University of California,Berkeley,Tech.Rep.UCB/EECS-2008-49,May2008” 和 “S.。

13、Tomov,R.Nath,H. Ltaief,and J.Dongarra,“Dense linear algebra solvers for multicore with gpu accelerators,” IEEE International Symposium on Parallel Distributed Processing Workshops and Phd Forum IPDPSW,pp.18,2010” 。 但是， 简单计算可知这样做的效率是极 低的。以稀疏矩阵 onetone2 （来自佛罗里达大学稀疏矩阵库） 为例， 该矩阵大小为 36k*36k， 对其进行稠密 LU 分解。

14、， 即使性能到达 1000Gflop/s， 仍需 15.5 秒， 而单核 CPU 上的稀疏 LU 分解只需要不到 1 秒。 0015 稀疏矩阵 LU 分解的方法的研究已经活跃了几十年。目前学术界与工业界已有许 多软件。稀疏矩阵 LU 分解的方法主要可分为两类。一类是基于稠密块的方法， 采用这类方 法的软件主要有 SuperLU 和 Pardiso （其中 Pardiso 有 GPU 版本） 。这类方法通过行列调换 在矩阵中形成稠密子块， 并对这些子块使用稠密矩阵的 LU 分解方法。但是， 在像电路矩阵 这样的极稀疏矩阵中， 通常无法形成稠密子块。因此， 这种方法不适合用于电路仿真。 0016 。

15、另一种是 G/P left-looking 方法， 该方法在 “J.R.Gilbert and T.Peierls,“Sparse partial pivoting in time proportional to arithmetic operations,” SIAM J.Sci.Statist.Comput.,vol.9,no.5,pp.862874,1988” 一 文 中 提出。该方法较强的数据依赖性使得它的并行化十分困难。据了解， G/P left-looking 方法尚未在 GPU 上实现。同时， 数据间的高依赖性还造成并行 G/Pleft-looking 方法只 能高效运行在共享内。

16、存的计算设备上， 如现场可编程门阵列 (Field Programmable Gate Array,FPGA)、 多核CPU和图形处理器(GPU)。 FPGA上实现的稀疏LU分解对大矩阵的可拓展 性受限于 FPGA 很有限的片上资源。由于共享内存的 CPU 核数通常很有限 （大多数商用 CPU 都不超过6核， 如Intel Xeon X5680,AMD Phenom II） ， 稀疏矩阵LU分解在多核CPU上性能 无法进一步提升。 发明内容 0017 本发明旨在至少解决上述技术问题之一。 说 明 书 CN 103399841 A 5 3/10 页 6 0018 为此， 本发明的目的在于提出一种。

17、基于 GPU 的稀疏矩阵 LU 分解方法， 该方法可根 据分解中浮点运算次数的多少自动选择计算平台， 能够合理分配CPU和GPU上的任务， 保证 GPU 上众多线程间的协同运行， 从而提高运算速度， 减少电路仿真时间。 0019 为了实现上述目的， 本发明的实施例提出了一种基于 GPU 的稀疏矩阵 LU 分解方 法， 包括以下步骤 ： A ： 在 CPU 上建立待仿真的电路网单的电路矩阵， 并对所述电路矩阵进行 预处理和稀疏化 ； B ： 根据预处理结果选择对所述稀疏矩阵 LU 进行处理的计算平台 ； C ： 如 果对所述稀疏矩阵 LU 进行处理所使用的计算平台为 GPU， 则通过所述 GPU。

18、 对所述稀疏矩阵 LU进行分解， 包括 ： C1 ： 将所述电路矩阵的非零元位置、 值以及所述稀疏矩阵的LU因子的非 零元位置输入所述 GPU ； C2 ： 确定群并行模式和流水线并行模式的分界点 ； C3 ： 依次计算属 于所述群并行模式的各层直至所有层均计算完成 ； C4 ： 通过所述流水线并行方式计算剩余 各层 ； C5 ： 将所述稀疏矩阵LU中非零元的值写回到所述CPU ； D ： 通过所述CPU对所述稀疏矩 阵 LU 中非零元的值进行回代求解， 以完成对稀疏矩阵 LU 的分解。 0020 根据本发明实施例的基于GPU的稀疏矩阵LU分解方法， 可根据分解中浮点运算次 数的多少自动选择计。

19、算平台， 能够合理分配 CPU 和 GPU 上的任务， 保证 GPU 上众多线程间的 协同运行， 从而提高运算速度， 减少电路仿真时间 ； 另外， 对于分解中浮点运算次数较多的 矩阵， GPU 稀疏矩阵 LU 分解相对单核 CPU、 4 核 CPU 和 8 核 CPU 分别可实现 7.9 倍、 2.6 倍和 1.5 倍的加速比。 0021 另外， 根据本发明上述实施例的基于GPU的稀疏矩阵LU分解方法还可以具有如下 附加的技术特征 ： 0022 在本发明的实施例中所述对所述电路矩阵进行预处理进一步包括 ： 对数值进行处 理以提高所述数值的稳定性 ； 根据近似最小自由度算法减少分解过程中的非零元。

20、写入 ； 根 据部分选主元的数值分解计算所述矩阵 LU 的非零元结构。 0023 在本发明的实施例中， 根据 HSL_MC64 方法对所述数值进行处理以提高所述数值 的稳定性。 0024 在本发明的实施例中， 所述根据预处理结果选择对所述稀疏矩阵 LU 进行处理的 计算平台， 进一步包括 ： 对所述矩阵 LU 中每行的非零元按列号进行排序 ； 建立消去图， 并对 所述消去图进行分层 ； 根据浮点运算次数确定要使用的平台 ； 如果所述浮点运算次数大于 或等于预设次数， 则使用 GPU 平台 ； 如果所述浮点运算次数小于所述预设次数， 则使用 CPU 平台。 0025 在本发明的实施例中， 所述建。

21、立消去图， 并对所述消去图进行分层， 进一步包括 ： 建立消去图， 其中， 所述消去图中包括节点， 且各节点表示各行。 对所述消去图进行分层， 其 中， 所述节点的层号为所述节点所依赖的行所处的最大层号加 1。 0026 在本发明的实施例中， 所述对所述消去图进行分层之后， 还包括 ： 依次处理所述消 去图的每一层， 所述每一层内部的各行可按照任意顺序进行计算。 其中， 处理所述消去图的 方式包括 ： 群并行模式和流水线并行模式， 所述群并行模式指所述每一层有足够多的行被 并行分解 ； 所述流水线并行模式指所述每一层只有很少的行， 且来自不同层的行重叠执行 分解。 0027 在本发明的实施例中。

22、， 所述通过所述流水线并行方式计算剩余各层， 进一步包括 ： 如果检测到未计算完成的行， 则设为第 k 行， 并计算所述第 k 行 ； 依次检测所述第 k 行所依 说 明 书 CN 103399841 A 6 4/10 页 7 赖的行， 并设为第j行， 并判断所述第j行是否计算完成 ； 如果所述第j行计算完成， 则从第 k 行中减去第 j 行乘以其相应系数 ； 如果第 j 行未完成， 则等待预设时间后重新检测第 j 行 是否计算完成 ； 依次确认所述第k行所依赖的每一行均计算完成后， 从所述第k行中减去该 行乘以其相应系数 ； 将第 k 行的计算结果写回到所述 GPU 中， 并判定第 k 行已。

23、计算完成， 并 进一步处理下一个未完成的行， 直至所有行均已完成计算。 0028 本发明的附加方面和优点将在下面的描述中部分给出， 部分将从下面的描述中变 得明显， 或通过本发明的实践了解到。 附图说明 0029 本发明的上述和 / 或附加的方面和优点从结合下面附图对实施例的描述中将变 得明显和容易理解， 其中 ： 0030 图 1 为根据本发明一个实施例的基于 GPU 的稀疏矩阵 LU 分解方法的流程图 ； 0031 图 2 为根据本发明一个实施例的基于 GPU 的稀疏矩阵 LU 分解方法的一个消去图 示意图 ； 0032 图 3 为根据本发明另一个实施例的基于 GPU 的稀疏矩阵 LU 分。

24、解方法的流程图 ； 0033 图 4 为根据本发明一个实施例的基于 GPU 的稀疏矩阵 LU 分解方法的对非零元访 问情况示意图 ； 0034 图 5 为根据本发明一个实施例的基于 GPU 的稀疏矩阵 LU 分解方法的不活跃 warp 引起死锁的示意图 ； 0035 图 6 为根据本发明一个实施例的基于 GPU 的稀疏矩阵 LU 分解方法的基于 GTX580GPU 的稀疏矩阵 LU 分解实现的带宽与活跃 warp 数之间的关系示意图 ； 和 0036 图 7 为根据本发明一个实施例的基于 GPU 的稀疏矩阵 LU 分解方法的各矩阵的浮 点运算次数及到达的带宽示意图。 具体实施方式 0037 下。

25、面详细描述本发明的实施例， 所述实施例的示例在附图中示出， 其中自始至终 相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。 下面通过参考附 图描述的实施例是示例性的， 仅用于解释本发明， 而不能理解为对本发明的限制。 0038 在本发明的描述中， 需要理解的是， 术语 “中心” 、“纵向” 、“横向” 、“上” 、“下” 、“前” 、 “后” 、“左” 、“右” 、“竖直” 、“水平” 、“顶” 、“底” 、“内” 、“外” 等指示的方位或位置关系为基于 附图所示的方位或位置关系， 仅是为了便于描述本发明和简化描述， 而不是指示或暗示所 指的装置或元件必须具有特定的方位、 以特。

26、定的方位构造和操作， 因此不能理解为对本发 明的限制。此外， 术语 “第一” 、“第二” 仅用于描述目的， 而不能理解为指示或暗示相对重要 性。 0039 在本发明的描述中， 需要说明的是， 除非另有明确的规定和限定， 术语 “安装” 、“相 连” 、“连接” 应做广义理解， 例如， 可以是固定连接， 也可以是可拆卸连接， 或一体地连接 ； 可 以是机械连接， 也可以是电连接 ； 可以是直接相连， 也可以通过中间媒介间接相连， 可以是 两个元件内部的连通。对于本领域的普通技术人员而言， 可以具体情况理解上述术语在本 发明中的具体含义。 说 明 书 CN 103399841 A 7 5/10 页。

27、 8 0040 本发明描述中的 “warp” 指的是 GPU 上的线程组， 一个线程组中包含的线程数量可 以通过 GPU 的产品手册查到。在本发明的一个实施例中， 使用的 NVIDIA GTX580GPU 上一个 warp 包含 32 个线程。在本发明的描述中，“warp” 与 “线程组” 是等价的术语。 0041 以下结合附图详细描述根据本发明实施例的基于 GPU 的稀疏矩阵 LU 分解方法。 0042 图 1 为根据本发明一个实施例的基于 GPU 的稀疏矩阵 LU 分解方法的流程图。如 图 1 所示， 根据本发明一个实施例的基于 GPU 的稀疏矩阵 LU 分解方法， 包括以下步骤 ： 00。

28、43 步骤 S101， 在 CPU 上建立待仿真的电路网单的电路矩阵， 并对电路矩阵进行预处 理和稀疏化。 0044 具体地， 对电路矩阵进行预处理时， 首先对数值进行处理以提高数值的稳定性， 在 本发明一个实施例中， 通过 HSL_MC64 方法对数值进行处理以提高数值的稳定性， 然后根据 近似最小自由度算法减少分解过程中的非零元写入， 最后根据部分选主元的数值分解计算 矩阵 LU 的非零元结构。 0045 步骤 S102， 根据预处理结果选择对稀疏矩阵 LU 进行处理的计算平台。 0046 具体地， 对稀疏矩阵 LU 中每行的非零元按列号进行排序， 建立消去图， 并对消去 图进行分层， 然。

29、后根据浮点运算次数确定要使用的平台， 具体而言， 如果浮点运算次数大于 或等于预设次数， 则使用 GPU 平台， 如果浮点运算次数小于预设次数， 则使用 CPU 平台。其 中， 预设次数根据具体情形预先设定。 0047 另外， 在上述步骤 S102 中， 建立消去图， 并对消去图进行分层进一步包括 ： 首先建 立消去图， 其中， 该消去图中包括节点， 且各节点表示各行， 其次对消去图进行分层， 且各节 点的层号为该节点所依赖的行所处的最大层号加 1。 0048 进一步地， 在对消去图进行分层之后， 依次处理消去图中的每一层， 且每一层内部 的各行可按照任意顺序进行计算， 其中， 处理消去图的方。

30、式包括 ： 群并行模式和流水线并行 模式。群并行模式指每一层有足够多的行被并行分解 ； 流水线并行模式指每一层只有很少 的行， 且来自不同层的行重叠执行分解。 0049 作为具体的示例， 在上述步骤 S102 中， 换言之， 即首先对矩阵 L 和 U 中每行的非零 元按照列号排序， 以提高数据局域性。在 GPU 程序中， 如果连续的线程访问连续的内存空 间， 可以实现较高的带宽。但预处理后， LU 因子中的非零元是乱序的， 这影响了访存的连续 性。因此对 L 和 U 中每行的非零元按列号排序， 以提高主存中的数据的局域性。如图 4 所 示， 排序后， 一行中相邻的非零元更有可能被相邻的线程处理。

31、。从而， 无需对 GPU 程序做任 何修改， 就获得了 1.7 倍的加速。排序过程只需进行一次， 并将其融合到预处理中。实验证 明， 对非零元进行快速排序的代价很小。 0050 其次， 建立消去图， 并对消去图分层。消去图中各节点表示行， 当且仅 lik 0 时， 即第 i 行依赖于第 k 行， 也即需要从 xi 中减去 uk， 生成从节点 k 指向节点 i 的边。如图 2 所示， 为一个消去图的示意图。其中， 任何符合消去图的更新顺序都是可行。具体来说， 只 要保证对任意的两行 xi 和 uk（xi 依赖于 uk） ， 从 xi 中减去 uk 时， uk 已经计算完成， 这个 顺序就是可行的。

32、。如图 2 所示， 是一个消去图的例子， 消去图可表示 left-looking 方法中 的并行性和要求的时序关系。假设图中虚线圆圈表示的行已经分解完成， 而行 8,9,10 正在 被处理。实线箭头代表的操作是目前可以执行的， 这些操作中是存在并行性的。同时要注 意， 对某一未完成行的更新操作需要按照严格的顺序执行。例如， 在图 2 所示的时刻， 行 9 说 明 书 CN 103399841 A 8 6/10 页 9 可以依次减去行 4,6,7， 对应图 2 中的实线箭头。而且， 行 9 减去行 4,6,7 后， 还必须等待行 8 完成， 才能减去行 8。对行 10 也存在类似的情况。进一步地。

33、， 对消去图分层， 某个节点的 层号， 是它所依赖的行所处的最大层号加 1。依次处理每一层， 每层内部的各行可以按照任 意顺序计算。 这种顺序是可行的， 因为某一行所依赖的行必定都处于它之前的层， 都会先于 该行被计算。消去图中， 属于同一层的行可以被独立地并行地处理。然而， 除了最前面的几 层以外， 后面每层往往只含有一两个节点， 不具备足够的并行性。此时， 需要将不同层的计 算重叠起来。 为此， 将每层有足够多的行被并行的分解称为群并行模式， 将每层只有很少的 行， 来自不同层的行重叠执行称作流水线并行模式。 0051 最后， 根据分解过程中涉及的浮点运算次数确定使用的平台。 具体地， 如。

34、果浮点运 算次数大于或等于预设次数 F0， 则使用 GPU 平台 ； 反之， 使用 CPU 平台。其中， 预设次数 F0 的具体值与所使用的 CPU 和 GPU 型号有关。在本发明一个优选实施例中， 使用的是 Intel E5680CPU 和 NVIDIA GTX580GPU， 相应的 F0=200M。 0052 步骤 S103， 如果对稀疏矩阵 LU 进行处理所使用的计算平台为 CPU， 则采用中国发 明专利 “201110337789.6， 陈晓明， 汪玉， 杨华中， 针对电路仿真的自适应并行 LU 分解方法” 中的方法进行稀疏 LU 分解 ； 如果对稀疏矩阵 LU 进行处理所使用的计算平。

35、台为 GPU， 则通过 GPU 对稀疏矩阵 LU 进行分解， 具体包括以下步骤 : 0053 步骤 S1031， 将电路矩阵的非零元位置、 值以及稀疏矩阵的 LU 因子的非零元位置 输入 GPU。 0054 步骤 S1032， 确定群并行模式和流水线并行模式的分界点。具体地， 从消去树的第 一层开始逐层检查每层上的点的个数， 直到某一层上的点数小于预设的 warp（线程组） 数 量时， 这一层之前的所有层采用群并行模式， 这一层以及这一层之后的所有层采用流水线 模式， 所述分界点确定。 0055 步骤 S1033， 依次计算属于群并行模式的各层直至所有层均计算完成。具体地， 依 次计算属于群并。

36、行模式的各层， GPU 上的每个线程组负责计算一行， 一个工作组内部的线程 操作同一行的不同非零元。具体过程即， 从正在被计算的行 （设为第 k 行） 中依次减去它所 依赖的行 （设为第 j 行） 乘以其相应系数 （即 lkj） 。通过维护一个标记数组， 记录每行是否 已完成计算， 初始化为未完成， 某行计算结束后， 标记该行为已完成， 最后重复上述过程直 到所有属于群并行模式的层计算完成。 0056 步骤 S1034， 通过流水线并行方式计算剩余各层。具体地， 如果检测到未计算完成 的行， 则将该行设为第 k 行， 并计算第 k 行， 依次检测第 k 行所依赖的行， 并设为第 j 行， 并 。

37、判断第 j 行是否计算完成， 如果第 j 行计算完成， 则从第 k 行中减去第 j 行乘以其相应系数 （ljk） ， 如果第j行未计算完成， 则等待预设时间后重新检测第j行是否计算完成 （此时一定 有某个线程组正在并行地计算第 j 行， 第 j 行最终一点会被完成， 因而此处不会形成死锁） ， 重复上述过程， 即依次确认第k行所依赖的每一行均计算完成后， 则从第k行中减去该行乘 以其相应系数， 最后将第 k 行的计算结果写回到 GPU 中， 并判定第 k 行已计算完成， 并继续 处理下一个未完成的行， 直至所有行均已完成计算， 即矩阵的 LU 分解已完成。 0057 步骤 S1035， 将稀疏。

38、矩阵 LU 中非零元的值写回到 CPU。即将 GPU 中的分解结果写 回到 CPU。 0058 步骤 S104， 通过 CPU 对稀疏矩阵 LU 中非零元的值进行回代求解， 已完成对稀疏矩 说 明 书 CN 103399841 A 9 7/10 页 10 阵 LU 的分解。 0059 在上述的示例中， 涉及到了GPU众多线程之间的协同， 以及保证GPU不同线程组之 间的时序关系， 这是基于 GPU 的稀疏矩阵 LU 分解方法中的一个难点。在 GPU 上保证时序关 系必须精确控制发布的线程组的个数。普通 GPU 程序中允许一些 warp 在开始时不活跃， 而 等其他 warp 完成后再活跃。然而。

39、， 在并行 left-looking 算法的流水线并行模式中， 我们必 须确保发布的所有 warp 从一开始就活跃， 否则将会造成死锁。如图 5 所示， 说明了这一情 况。假设在一个只支持两个活跃 warp 的 GPU 上发布了三个 warp， 群并行模式可以顺利完 成， 因为 warp1 和 2 最终总会结束运行， 于是 warp3 开始运行。但在流水线并行模式中， 行 9 和行 10 依赖于行 8， 而行 8 被分配给了不活跃的 warp3， 于是活跃的 warp（1 和 2） 一直等 待， 陷入死循环。反过来又使 warp3 没有机会活跃， 从而形成死锁。因此流水线并行模式中 能被并行分。

40、解的最大行数就是 GPU 程序执行时活跃 warp 的数目。 0060 对许多 GPU 程序而言， 更多的活跃 warp 数目通常意味着更高的性能， 但如果发布 了太多的活跃 warp， 同时被分解的行数过多， 可能会出现如下情况 ： 矩阵下方的行依赖于 一些之前的行， 由于这些行没有完成， 下方的这些行只能等待。 在这里等待使用死循环实现 的， 这些行的分解无法进行， 却仍然会占用运算单元。这反而会造成性能的下降。这个转折 点不会在 CPU 上出现， 但在 GTX580GPU 上可以观察到这个现象。如图 6 所示， 在图 6 中， 给 出了在 GTX580 总分解其中四个矩阵， 在发布不同 。

41、warp 数时达到的带宽。最高性能并不是 在活跃 warp 数最多时达到的， 而是在大约每个流处理器有 24 个 warp 时达到的。这意味着 在 GTX580 中已经充分利用了稀疏矩阵 LU 分解中流水线并行模式中的并行性。采用一个自 动调整策略来确定 GPU 上最佳的 warp 数， 即针对某个输入矩阵， 使用不同的 warp 数目， 多 次执行上述 GPU 稀疏矩阵 LU 分解方法， 找到使性能最佳的 warp 数， 用于后续迭代之中。 0061 在 GTX580 上， 最高可实现 74% 的峰值带宽。考虑到访存并不完全连续， 而且 warp 有时候需要等待， 所以这个带宽已经接近峰值的。

42、192GB/s。 如果今后GPU峰值带宽能够进一 步提升， 则稀疏矩阵 LU 分解的性能也将进一步提升。 0062 图 3 为根据本发明另一个实施例的基于 GPU 的稀疏矩阵 LU 分解方法的流程图。 0063 如图 3 所示， 根据本发明另一个实施例的基于 GPU 的稀疏矩阵 LU 分解方法， 包括 以下步骤 ： 0064 步骤 S301， 在 CPU 上进行预处理。其中， 预处理包括 ： 1） 使用 HSL_MC64 方法提高 数值稳定性 ； 2） 使用近似最小自由度算法减少分解过程中的非零元填入 ； 3） 使用部分选主 元的数值分解计算 L 和 U 的非零元结构 ； 4） 对矩阵 L 和。

43、 U 中每行的非零元按列号排序以提 高数据局域性。预处理在 CPU 上进行， 后续的迭代的数值分解过程在 GPU 上进行。在电路 仿真的迭代过程中， L 和 U 的非零元结构保持不变。将预处理后的矩阵记为矩阵 A。 0065 步骤 S302， 对矩阵 L 和 U 中的非零元进行排序。具体地， 对矩阵 L 和 U 中每行的非 零元按列号排序， 可提高数据局域性。 0066 进一步地， 通过对矩阵进行图的深度优先遍历， 得到消去图， 并将消去图分层。并 根据分解过程中涉及的浮点次数确定是否使用 GPU 平台。具体地， 如果浮点运算次数少于 F0， 则使用 CPU 平台 ； 如果浮点运算次数多于 F。

44、0， 则使用 GPU 进行稀疏矩阵 LU 分解， 并进一 步执行步骤 S303。 0067 步骤 S303， 将矩阵 A 的 L 和 U 中的非零元位置写入 GPU 内存。 说 明 书 CN 103399841 A 10 8/10 页 11 0068 步骤 S304， 将矩阵 A 的非零元的值写入 GPU 内存。 0069 步骤 S305， 判断是否处于群并行模式， 如果是则执行步骤 S306， 否则执行步骤 S307。 0070 步骤 S306， 在 GPU 上并行地分解该层中的各行。完成后返回继续执行步骤 S305。 具体而言， 在群并行模式中， 每层有足够多的行， 可以独立地并行地被分解。

45、。依次并行地处 理每层中的各行， 各行被分配给 GPU 的各个工作组， 一个工作组内部的线程操作同一行的 不同非零元。维护一个标记数组， 记录每行是否已完成计算。初始化为未完成， 处理完某行 后， 标记该行已完成。 0071 步骤S307， 在GPU上按流水线并行模式分解剩余所有行。 具体地， 在流水线并行模 式中， 每层的行很少， 其并行模式与群并行模式有重要区别。在流水线并行模式中， 需要保 证时序关系。每次更新前， 必须检查被依赖的行的计算是否已经完成。如果尚未完成， 则必 须等待一段时间之后重新检查， 用这种方法可保证所要求的时序关系。 0072 步骤 S308， 从 GPU 中读取矩。

46、阵 LU 中非零元的值， 进一步地， 将该值写回到 CPU， 并 回代求解。 0073 需要说明的是， 在上述图3所示的流程图的步骤中， 步骤S301至步骤S303， 只需执 行一次， 而步骤 S304 至步骤 S308 需要多次循环执行。 0074 作为具体的示例， 下表 1 显示了根据本发明一个实施例的基于 GPU 的稀疏矩阵 LU 分解的性能， 以及与 CPU 版本实现的加速效果。表 1 中所列测试矩阵来自佛罗里达大学矩 阵收集 （http:/www.cise.ufl.edu/research/sparse/matrices/） 。 0075 表 1.GPULU 分解性能及加速比 1 0。

47、076 说 明 书 CN 103399841 A 11 9/10 页 12 0077 1 除 KLU 之外， 其他 CPU 结果均指我们自己的 CPU 实现 0078 2 矩阵大小 0079 3 矩阵 A 中的非零元 0080 4 分解中的浮点运算次数 （兆为单位） 0081 5 所有的平均值均为几何平均 0082 表 1 0083 图 7 为根据本发明一个实施例的基于 GPU 的稀疏矩阵 LU 分解方法的各矩阵的浮 点运算次数及到达的带宽示意图。 0084 如图 7 所示， 给出了各个矩阵分解中的浮点运算次数以及分解该矩阵时在 GPU 上 实现的带宽。基于 GPU 的稀疏矩阵 LU 分解， 。

48、对于分解过程中浮点运算次数较多 （本发明的 实施例中大于 200M， 即 B 组） 的矩阵非常有效。在这些矩阵上， 基于 GPU 的 LU 分解相对单 核 CPU、 4 核 CPU、 8 核 CPU 分解具有 7.9 倍、 2.6 倍和 1.5 倍的加速。由于大多数商用 CPU 都 不超过 6 核， 因此对于分解中运算量较大的矩阵使用 GPU 进行 LU 分解， 能带来可观的性能 提升。且应用在电路仿真领域时， 还能提高 SPICE 电路仿真器的性能。 0085 需要说明的是， 在分解C组矩阵时， 性能大大超过CPU。 这是因为， 分解C组矩阵的 过程中会出现许多反常浮点数 （denormal。

49、 floating point numbers） 。 其中， 反常浮点数用 来表示极小的实数。而 CPU 处理这种数据很慢， 因此分解这些矩阵时性能很差。 0086 另外， GPU 上的稀疏矩阵 LU 分解本身也具有重要意义。目前， GPU 上的基本线性代 说 明 书 CN 103399841 A 12 10/10 页 13 数程序库(Basic Linear Algebra Subprogram， BLAS)中已有三角矩阵求解的模块， 但尚没 有稀疏矩阵 LU 分解模块。而 LU 分解的作用可将一个一般矩阵方程转换为两个三角矩阵方 程。 0087 根据本发明实施例的基于GPU的稀疏矩阵LU分解方法， 可根据分解中浮点运算次 数的多少自动选择计算平台， 能够合理分配 CPU 和 GPU 上的任务， 保证 GPU 上众多线程间的 协同运行， 从而提高运算速度， 减少电路仿真时间 ； 另外， 对于分解中浮点运算次数较多的 矩。