基于 AMD GPU 的稀疏矩阵 LU 分解 及其在电路仿真中的应用.pdf

基于基于 amd gpu 的稀疏矩阵的稀疏矩阵 lu 分解分解 及及其在电路仿真中的应用其在电路仿真中的应用 任令任令 陈晓明陈晓明 张晨曦张晨曦 指导教师指导教师 汪玉汪玉 摘要 稀疏矩阵求解已成为spice仿真器的瓶颈 由于稀疏矩阵lu分解中数据间的依赖性强 目前关于稀 疏矩阵求解的工作较少 较强的数据依赖性也决定了并行稀疏矩阵lu分解算法只有在共享内存的计算设备 上才能高效运行 但是共享统一内存的cpu核的个数通常很有限 最新的图形处理器 gpu 提供了一个 可能的解决方案 gpu拥有共享内存的众核结构 本文中 我们首次提出了一个基于gpu的稀疏矩阵lu求 解器 我们针对gpu架构优化了任务划分 活跃线程数和内存访问模式 对于分解中浮点运算次数较多的 矩阵 我们的gpu稀疏矩阵lu分解相对单核cpu和8核cpu分别实现了10 5x和2 6x的加速比 一 一 引言引言 spice simulation program with integrated circuit emphasis 1 是当今最常用的电路仿真器 然而 超大规模集成电路 vlsi 的快速发展对spice的性能提出了挑战 在当今的超大规模集 成电路中 后仿真寄生参数提取后的电路矩阵的尺寸很容易到达几兆 这使得后仿真在cpu上需 要几天甚至数周的时间才能完成 spice仿真中两个最耗时的步骤是利用lu分解求解稀疏矩阵 方程和参数更新 这两个步骤要被迭代进行许多次 近十年来 gpu在许多领域得到了重要应用 2 3 目前已有一些工作使用gpu完成参数更 新 4 5 参数更新的并行化非常直观 它由一系列独立的任务组成 即电路中各个元件的参数 更新 然而 由于稀疏矩阵的结构很不规则 且lu分解中的数据间的依赖性很强 稀疏矩阵求 解器的并行化十分困难 据我们了解 目前尚没有能适合电路仿真的基于gpu的稀疏矩阵求解器 以往的gpu lu分解工作集中在稠密矩阵上 6 9 尽管基于gpu的稠密矩阵lu分解性能非 常可观 在gtx 280上达到388 gflop s giga floating point operations per second 8 简单计 算可知我们不应该将稀疏矩阵当作稠密矩阵分解 以onetone2 大小为36k 36k 为例 即 使将稠密lu分解的性能计算为1000 gflop s 对onetone2进行稠密lu分解仍需要15 5s 而单 核cpu上一个简单的串行稀疏lu分解只需要不到1s 目前 稀疏矩阵lu求解器主要在cpu上实现 为了利用更多的并行性 superlu 10 12 将 supernode思想与gilbert peierls g p 的leftlooking 13 算法结合起来 实现了对具有稠密块的矩 阵的高效分解 pardiso 14 也利用了supernode christen等人在gpu上实现了pardiso 15 他们的思路仍然遵循基于gpu的稠密矩阵分解 他们用gpu计算稀疏矩阵中的稠密块 而其余的 工作在cpu上完成 然而 supernode在电路矩阵等极稀疏的矩阵中很少出现 这使得基于supernode的算法在处 理电路矩阵时比基于列的算法效率低 所以在klu这款针对电路仿真优化的稀疏求解器中 采 用了直接的g p left looking算法 13 和对角块矩阵形式block triangular form btf 16 而放弃 了supernode思想 klu只有串行版本 由于数据间高度依赖 并行g p leftlooking算法只能高效运行在共享内 存的计算设备上 如现场可编程门阵列 field programmable gate array fpga 多核cpu and gpu 一些工作在fpga实现了g p leftlooking算法 17 19 但这些方法对大矩阵的可拓展性还 受限于fpga很有限的片上资源 我们在多核cpu上也实现了这一算法 20 但能共享内存的cpu 核数通常很有限 大多数商用cpu都不超过6核 如intel xeon x5680 amd phenom ii 最新的 gpu为该问题提供了一个可能的解决途径 与cpu相比 gpu共享内存的核数远远多于cpu 且 拥有更高的内存带宽 因此在本文中 我们首次针对电路仿真提出了一个基于的gpu的不选主元 的稀疏矩阵lu分解器 我们通过以下优化更好的利用了算法中的并行性 挖掘更多挖掘更多并行性并行性 gpu拥有成千上万的同时运行的线程 所以我们必须挖掘足够多的 并行性从而在gpu实现高效的稀疏矩阵lu分解器 在 20 中 我们提出了两种并行模 式 来描述向量操作之前的并行性 但是仅这一层并行性对于gpu来说是不够的 我 们进一步利用了向量操作内部的并行性 为了有效的处理这两层并行性 我们根据两 种模式的不同特点和gpu的架构划分任务 我们的策略减少了空闲线程 节省了线程 同步的开销 并保证了全局有足够的线程数 保证时序关系 保证时序关系 并行left looking算法中必须保证合适的时序关系 在gpu上保证时序 关系涉及到精确控制线程组 amd gpu中的wavefront 的数目 优化访存模式 优化访存模式 我们 1 设计了适合gpu的中间结果数据格式 2 提出对非零元 排序以提高数据局域性 使得对global memory的访问更连续 在36个矩阵上的实验结果说明基于gpu的稀疏矩阵lu分解 对于分解过程中浮点运算次数 较多 我们的平台中大于600m 的矩阵非常有效 在这些矩阵上 gpu相对单核cpu和8核cpu 分别实现了10 5x和2 6x的加速比 对于分解中出现反常数 21 的矩阵加速比更大 本文其余部分结构如下 第二章详细介绍我们的基于gpu的稀疏矩阵lu分解 第三章给出 实验结果和讨论 第四章结束全文 二 二 基于基于 gpu 的稀疏矩阵的稀疏矩阵 lu 分解分解 图1 基于gpu的稀疏矩阵lu分解流程图 图1是我们基于gpu的稀疏矩阵lu分解的流程图 预处理只在cpu上进行一次 第二章第1 节 数值分解部分 我们挖掘g p left looking算法中更多的并行性 并根据两种并行模式的不 同特点和gpu的架构划分任务 第二章第2节 然后我们讨论基于gpu的稀疏矩阵lu分解算法 中的几个关键问题 包括gpu线程组之间的时序关系 第二章第3节 和对访存模式的优化 第 二章第4节 1 预处理预处理 预处理包括3个步骤 1 为了提高数值稳定性 使用hsl mc64算法 22 2 为了减少分 解过程中的非零元填入 使用近似最小自由度算法 amd approximate minimum degree 23 3 为了计算l和u的非零元结构 使用基于g p算法的符号分解 部分选主元的数值分解 13 在 电路仿真的迭代过程中 l和u的非零元结构保持不变 我们记预处理后的矩阵为a 2 并行性挖掘并行性挖掘 算法1是串行的g p left looking算法 其中的核心操作时向量乘加 mad 我们在 20 提出了 该算法中的两种并行模式 cluster mode 和 pipeline mode 这两个并行模式描述了向量乘加操作 之间的并行性 我们以图2为例说明这一层并行性 假设途中虚线圆圈表示的列已经分解完成 而列8 9 10正在被处理 用绿色实线箭头表示的操作都是可以执行的 虽然对对某一未完成列的 更新操作需要按照严格的顺序执行 但这些操作中仍然是存在并行性的 但是仅靠向量乘加操作 间的并行性还不能利用gpu的高带宽 我们利用了算法中另一层内在的并行性 向量操作内部的 并行性 为了有效的处理这两层并行性 我们通过任务划分充分利用gpu资源 这个过程中 有 几个关键因素需要考虑 为了方便 无论我们用多少个线程来分解一列 我们把它们称作一个虚拟组 同一个虚拟 组中的线程操作同一列中的非零元 且必行同步 较小的虚拟组可以减少空闲线程 但是基于以 下两个原因 虚拟组也不应太小 首先 同时能被分解的列数 即虚拟组个数受到存储空间的限 制 所以太小的虚拟组会导致全局的线程数过少 这对gpu程序中不希望的 第二 gpu按照单 指令多数据的方式 simd single instruction multiple data 调度线程 一组simd线程在amd gpu 上被称为一个wavefront 如果同一个wavefront内部的线程出现分支 所必要的路径将被串行地 执行 不同的虚拟组处理不同的列 从而经常出现分支 因此从减少simd线程分支的角度考虑 一个虚拟组至少应该包括一个wavefront 图2 并行left looking算法中的并行性和时序关系 综合考虑以上因素 我们提出下列任务划分策略 cluster mode中的列非常稀疏 因此我们 在保证全局线程数足够多前提下 尽量减小虚拟组的大小 而pipeline mode中的列通常非零元较 多 足够一个或几个wavefronts来处理 此时虚拟组的大小在减少空闲线程的意义上并不重要 我们选择用一个wavefront处理一列 这样不仅减小了simd线程之间的分支 也减少了线程同步 的开销 因为一个wavefront的线程间的同步已经自动由gpu的simd架构保证了 3 gpu 上时序关系的保证上时序关系的保证 算法2是pipeline mode的并行left looking算法 这种并行模式必须保证列之间恰当的时序关系 如果列k依赖于列t 只有当列t完成后 列k才能被列t更新 我们仍然用图2的例子来解释所需的 时序关系 假设列8 9 10正在处理 其他列已经完成 此时列9可以首先被列4 6 7 更新 对应 途中的绿色实线箭头 但目前列9不能被列8更新 它必须等待列8完成 对列10也存在类似的情 况 在gpu上保证时序关系值得特别注意 此时 gpu程序中wavefront的个数必须精确控制 这与gpu活跃wavefront的概念有关 一个gpu kernel可以有很多的wavefront 通常 由于资 源限制 一些wavefront不能立即开始执行 而要等待其他wavefront完成后才能变为 活跃 然而 在并行稀疏矩阵lu分解的pipeline mode中 我们必须确保我们发布的所有wavefront 从一开始就是活跃的 设想某一列被分配给了一个不活跃的wavefront 依赖于这一列的 列就必须等待它分解完成 但反过来 负责这个未完成的列的wavefront一直没有机会 活跃 因为当前活跃的wavefront都无法结束运行 这将形成死锁 图3说明了这一情况 假设我们在一 个只支持两个活跃wavefront的gpu上发布了三个wavefront cluster mode中不会出现问题 因为 wavefront 1和2最终总会结束运行于是wavefront 3可以开始 但在pipeline mode中 列9和列10依 赖于列8 而列8被分配给了一个不活跃的wavefront 于是活跃的wavefront 1和2 陷入死循环 一直等待列8 这反过来又是的wavefront 3没有机会活跃 从而形成死锁 图3 不活跃wavefront引起死锁的说明 因此 能被pipeline mode中并行分解的最大列数恰好等于kernel执行时活跃wavefront的数目 这个数目依赖于许多因素 如资源使用量和分支与循环数目 24 并显著影响我们的gpu稀疏lu 分解器的性能 4 访存模式的优化访存模式的优化 gpu稀疏矩阵lu分解的优化主要是针对访存模式的优化 本小节讨论 1 我们对中间结果 向量数据格式的选择 和 2 对非零元排序以提高数据局域性 使得对global memory的访问更 连续 中间向量格式 中间向量格式 我们有两种可选的方式来储存中间结果 即算法2中的x csc compressed sparse column 稀疏向量或稠密数组 csc稀疏向量节省空间 可以存储在local memory中 而 稠密数组只能存储在global memory中 通常 尽可能地利用local memory可以提高性能 但在本 问题中却不是这样 首先csc格式的下标访问很不方便 要访问csc格式的稀疏向量中的元素 最快的方法就是二分搜索 binary search 即使访问的是local memory 效率仍然很低 更重 要的是 过多使用local memory会减少每个运算单元 compute unit cu 上的活跃wavefront的 数目 每个cu的活跃work group数 每个cu的local memory大小 一个csc稀疏向量的大小 从而使性能严重下降 因此我们采用稠密数组储存中间结果 提高数据局域性 提高数据局域性 在gpu中 如果连续的work item访问连续的内存空间 可以实现较高的 带宽 但预处理后 lu因子中非零元是乱序的 这影响了访存的连续性 我们用快速排序算法 对l和u中的按照行号排序 以提高主存中的数据的局域性 如图4所示排序后 一列中相邻的非 零元更有可能被相邻的线程处理 排序过程只需进行一次 我们将它融合到预处理中 实验证明 对非零元排序的代价很小 而可以使原kernel加速1 7倍 图4 排序之后访问global memory更连续 三 三 实验结果实验结果与讨论与讨论 1 实验环境实验环境 cpu程序用c编写 运行环境为64位linux系统 两块xeon e5405 cpu 共8核 具有 2 6 mb l2缓存 16gb内存 我们的cpu实现利用了sse streamingimd extensions gpu程 序用opencl v1 1 25 编写 运行环境为amd app sdk v2 4 amd radeon 5870 gpu 我们用 36个来自佛罗里达大学稀疏矩阵库 26 的矩阵进行测试 虽然我们的应用背景是电路仿真 我们 也包括了一些来自其它应用的稀疏矩阵 以证明我们的稀疏矩阵lu分解器并不局限于电路矩阵 2 性能与加速比性能与加速比 表i给出了gpu稀疏矩阵lu分解的性能 并与不同核数的cpu性能进行比较 表中所列时间 是数值分解时间 不包括预处理和回代求解的时间 有些数据传输在每次分解中都要进行 这些 时间被计入了gpu运行时间中 我们发现gpu的加速比与分解某矩阵所需的浮点运算次数 flops 密切相关 所以 全部矩阵平均加速比和带宽并不反映太多有用信息 我们将所用的测试矩阵分 为四组 前三组根据浮点运算次数分组 表1 gpu lu分解性能及加速比1 matrix 2n k 3non zeros k 4mop 1 core cpu time s 1 core cpu bandwidth gb s gpu time s gpu bandwidth gb s speedup over 1 core 4 core 8 core klu hcircuit 103 2 513 1 1 0 0 02 1 38 0 04 0 70 0 51 0 32 0 32 0 70 lung2 109 5 492 6 1 1 0 01 2 24 0 03 0 91 0 41 0 23 0 23 0 68 circuit 4 80 2 307 6 2 5 0 02 3 07 0 07 0 96 0 31 0 15 0 13 0 43 rajat21 402 0 1893 4 3 3 0 09 0 97 0 58 0 16 0 16 0 16 0 15 0 24 bcircuit 67 3 375 6 5 1 0 05 3 03 0 04 3 49 1 15 0 34 0 34 1 31 5avreage 1 1 95 0 81 0 41 0 23 0 22 0 58 dc1 116 8 766 4 16 9 0 12 3 91 0 48 0 96 0 25 0 09 0 07 0 25 trans4 116 8 766 4 16 9 0 12 3 83 0 48 0 96 0 25 0 09 0 07 0 26 hvdc2 189 9 1347 3 19 2 0 15 3 50 0 10 5 25 1 50 0 84 0 53 6 17 onetone2 36 1 227 6 94 1 0 49 5 24 0 08 32 16 6 14 1 84 1 02 12 54 transient 178 9 961 8 107 8 0 60 4 94 0 36 8 19 1 66 0 50 0 30 1 68 ckt11752 dc 1 49 7 333 0 144 6 0 76 5 21 0 18 21 97 4 22 1 30 0 85 0 55 tsopf rs b300 c3 42 1 4413 5 211 1 1 18 4 90 0 20 28 86 5 89 2 92 2 53 1 25 epb3 84 6 463 6 267 2 1 42 5 13 0 21 34 79 6 78 1 92 1 09 34 44 raj1 263 7 1302 5 340 7 1 87 4 97 0 87 10 71 2 15 0 73 0 48 220 39 asic 680ks 682 7 2329 2 436 5 3 05 3 91 0 49 24 36 6 23 1 99 1 21 5 86 thermomech tc 102 2 711 6 449 4 2 40 5 12 0 37 33 21 6 49 2 00 1 14 5 70 asic 680k 682 9 3871 8 474 8 3 53 3 68 2 00 6 49 1 76 0 60 0 39 2 19 asic 100k 99 3 954 2 529 6 2 95 4 91 0 69 20 99 4 28 1 35 0 86 6 06 avreage 2 4 51 10 97 2 43 0 84 0 55 3 64 asic 100ks 99 2 578 9 663 0 3 71 4 88 0 34 53 32 10 92 3 44 2 21 15 38 rma10 48 6 2374 0 730 6 3 81 5 25 0 57 35 05 6 68 2 28 1 25 227 25 onetone1 36 1 341 1 799 8 3 48 6 29 0 41 53 34 8 48 2 84 1 83 61 75 thermomech dm 204 3 1423 1 898 8 4 84 5 08 0 71 34 61 6 81 2 15 1 32 5 96 venkat50 62 4 1717 8 1043 5 6 29 4 53 0 61 46 78 10 32 3 30 2 04 8 14 zhao1 33 9 166 5 1737 1 9 64 4 93 0 79 60 13 12 20 4 15 2 84 10 96 thermomech dk 204 3 2846 2 3637 9 22 22 4 48 2 05 48 52 10 84 3 78 2 56 8 72 crashbasis 160 0 1750 4 3933 0 22 46 4 79 2 12 50 73 10 59 3 79 2 56 9 42 g2 circuit 150 1 726 7 4780 0 26 51 4 93 2 10 62 24 12 62 4 59 3 52 22 20 twotone 120 8 1222 4 5245 0 28 06 5 11 2 10 68 29 13 36 4 76 4 01 11 19 sme3dc 42 9 3148 7 5291 9 29 10 4 97 2 60 55 65 11 19 3 98 3 20 10 04 xenon1 48 6 1181 1 10066 3 55 59 4 95 3 90 70 58 14 25 5 31 4 51 12 73 helm2d03 392 3 2741 9 13331 8 73 67 4 95 6 18 58 99 11 92 4 75 4 27 10 62 avreage 3 4 99 52 59 10 53 3 65 2 58 15 58 denormal 89 4 1156 2 2387 1 13 18 4 95 1 10 59 34 11 99 4 06 2 68 10 69 torso2 1033 5 116 0 651 4 6 91 2 58 0 47 37 90 14 70 4 52 2 42 13 81 majorbasis 160 0 1750 4 3933 0 46 22 2 33 2 12 50 73 21 80 6 97 3 97 20 65 asic 320k 321 8 2635 4 584 1 39 13 0 41 1 17 13 65 33 45 8 84 4 61 33 45 asic 320ks 321 7 1827 8 651 6 44 63 0 40 0 41 43 46 108 86 29 19 15 19 108 03 avreage 4 1 37 36 82 26 86 8 01 4 48 25 61 1 除klu之外 其他cpu结果均指我们自己的cpu实现 2 矩阵大小 3 矩阵a中的非零元 4 分解中的浮点运算次数 兆为单位 5 所有的平均值均为几何平均 1 分解第一组的矩阵需要的运算次数很少 在这些矩阵上 gpu不能获得加速比 甚至 性能比单核cpu更差 可以注意到 对于这些矩阵 cpu和gpu的带宽都非常低 这 说明程序其他的开销占用了大多数运行时间 2 第二组的矩阵分解运算量适中 此时cpu带宽已基本稳定 gpu在这些矩阵上的性能 通常略优于单核cpu但比多核cpu差 3 分解第三组中的矩阵需超过600兆浮点运算 gpu对于这些矩阵的分解加速比非常可观 相对于单核cpu和8核cpu分别有10 5倍和2 5倍加速比 gpu带宽可以到达50 70gb s 但与峰值154gb s仍有一定距离 在图5中 我们画出了前三组矩阵的浮点运算次数和gpu加速比之间的关系 从图中我们可 以看到gpu加速比正相关于浮点运算次数 除了极少数例外点 我们可以根据某个矩阵分解的浮 点运算次数来推断gpu性能能否超过8核cpu 这个结果也说明 至少在稀疏矩阵lu分解中 只 有在问题规模足够大的情况下 gpu超强的浮点计算能力和高带宽才能被很好的利用 图5 各矩阵的浮点运算次数及加速比 第四类矩阵比较特殊 分解这些矩阵的过程中会出现许多反常浮点数 denormal floating point numbers 反常浮点数被用来表示极小的实数 cpu处理这种数很慢 8 因此分解这些矩 阵时性能很差 目前较新的gpu可以全速处理反常数 这是amd gpu在科学计算应用中的一个 优势 3 局限性局限性 我们的稀疏矩阵lu分解的一个主要局限在于 无法处理l u i中非零元太多的矩阵 这 是因为radeon 5870 gpu的主存相对来说较小 1gb 但目前最先进的gpu 如amd radeon 6990 已经具备4gb内存 这对于大多数矩阵来说已经足够了 随着gpu的发展 我们的gpu稀疏矩阵 lu分解模块将能够处理具有更多非零元的矩阵 四 四 结论结论 本文首次提出了基于gpu的针对电路仿真的稀疏矩阵lu分解研究 我们详细分析了我们的 gpu稀疏矩阵lu分解的性能 并与cpu版本和klu做了比较 我们发现分解中浮点运算次数超 过较多的稀疏矩阵适合用gpu加速 对这些矩阵上 gpu相对单核cpu和8核cpu分别实现了 10 5x和2 6x的加速比 参考文献参考文献 1 l w nagel spice 2 a computer program to stimulate semiconductor circuits ph d dissertation university of california berkeley 1975 2 gpgpu 2 proceedings of 2nd workshop on general purpose processing on graphics processing units new york ny usa acm 2009 3 gpgpu 10 proceedings of the 3rd workshop on general purpose computation on graphics processing units new york ny usa acm 2010 4 n kapre and a dehon performance comparison of single precision spice model evaluation on fpga gpu cell and multi core processors in field programmable logic and applications international conference on 2009 pp 65 72 5 k gulati j f croix s p khatr and r shastry fast circuit simulation on graphics processing units in proceedings of the 2009 asia and south pacific design automation conference ieee press 2009 pp 403 408 6 n galoppo n k govindaraju m henson and d manocha lugpu efficient algorithms for solving dense linear systems on graphics hardware sc conference vol 0 p 3 2005 7 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 may 2008 8 s tomov j dongarra and m baboulin towards dense linear algebra for hybrid gpu accelerated manycore systems parallel comput vol 36 pp 232 240 june 2010 9 s 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 1 8 2010 10 j w demmel s c eisenstat j r gilbert x s li and j w h liu a supernodal approach to sparse partial pivoting siam j matrix analysis and applications vol 20 no 3 pp 720 755 1999 11 j w demmel j r gilbert and x s li an asynchronous parallel supernodal algorithm for sparse gaussian elimination siam j matrix analysis and applications vol 20 no 4 pp 915 952 1999 12 x s li and j w demmel superlu dist a scalable distributed memory sparse direct solver for unsymmetric linear systems acm trans mathematical software vol 29 no 2 pp 110 140 june 2003 13 j r gilbert and t peierls sparse partial pivoting in time proportional to arithmetic operations siam j sci statist comput vol 9 pp 862 874 1988 14 o schenk and k gartner solving unsymmetric sparse systems of linear equations with pardiso computational science iccs 2002 vol 2330 pp 355 363 2002 15 m christen o schenk and h burkhart general purpose sparse matrix building blocks using the nvidia cuda technology platform 2007 16 t a davis and e palamadai natarajan algorithm 907 klu a direct sparse solver for circuit simulation problems acm trans math softw vol 37 pp 36 1 36 17 september 2010 17 j johnson t chagnon p vachranukunkiet p nagvajara and c nwankpa sparse lu decomposition