大规模稀疏线性方程组的GMRES-GPU快速求解算法.pdf

第2 3 卷第4 期 2 0 1 1 年4 月 计算机辅助设计与图形学学报 J o u r n a lo fC o m p u t e r A i d e dD e s i g n & C o m p u t e rG r a p h i c s V 0 1 2 3N o 4 A p r 2 0 1 1 大规模稀疏线性方程组的G M R E S G P U 快速求解算法 柳有权1 ，尹康学”，吴恩华2 3 1 ( 长安大学信息工程学院西安7 1 0 0 6 4 ) 2 ( 中国科学院软件研究所计算机科学国家重点实验室北京1 0 0 1 9 0 ) 3 ( 澳门大学科技学院澳f - j ) ( y o u q u a n e h d e d u o n ) 擒耍：重开始广义极小残量法( G M R E S ) 是求解大规模线性方程组的常用算法之一，具有收敛速度快、稳定性好等 优点文中基于C U D A 将G M R E S 算法在G P U 上进行并行算法实现，尤其针对稀疏矩阵矢量乘法运算，通过合并访 问和共享内存策略相结合的手段使得算法效率大幅度提升对于大规模数据集，在G e F o r c eG T X2 6 0 上的运行结果 相对于I n t e lC o r e2Q u a dC P UQ 9 4 0 0 2 6 6 G H z 得到了平均4 0 余倍的加速效果，相对于I n t e lC o r ei 7C P U9 2 0 2 6 7G H z 也可得到平均2 0 余倍的加速效果 关键词：C U D A ；G P G P U ；重开始广义极小残量法；稀疏矩阵矢量乘法 中图法分类号：T P 3 9 1 F a s tG M R E S G P US o l v e rf o rL a r g eS c a l eS p a r s eL i n e a rS y s t e m s L i uY o u q u a n lt 幻，Y i nK a n g x u e ”，a n dW uE n h u a 2 3 ”( S c h o o lo fI n f o r m a t i o nE n g i n e e r i n g C h a n g a nU n i v e r s i t y ，X i a n7 1 0 0 6 4 ) ”( S t a t eK e yL a b o r a t o r yo fC o m p u t e rS c i e n c e 。I n s t i t u t eo fS o f t w a r e C h i n e s eA c a d e m yo fS c i e n c e s ，B e O i n g1 0 0 1 9 0 ) ”( F a c u l t yo fS c i e n c ea n dT e c h n o l o g y 。U n i v e r s i t yo fM a c a u ，M a c a o ) A b s t r a c t ：A sap o p u l a ri t e r a t i v em e t h o dt os o l v el i n e a re q u a t i o n s ，r e s t a r t e dg e n e r a l i z e dm i n i m a l r e s i d u a lm e t h o d ( G M R E S ) h a st h ea d v a n t a g e so ff a s tc o n v e r g e n c ea n dg o o ds t a b i l i t y T h i sp a p e r i m p l e m e n t sap a r a l l e lG M R E Si n G P Ub a s e dO i lC U D A P a r t i c u l a r l y ，t h es p a r s em a t r i xv e c t o r m u l t i p l i c a t i o ni so p t i m i z e dw i t hc o h e r e n c ev i s i t i n ga n ds h a r e dm e m o r y ，w h i c hs i g n i f i c a n t l yi m p r o v e s t h ep e r f o r m a n c e W et e s t e dt h ep a r a l l e l e dG M R E So naG P Uo fG e F o r c eG T X 2 6 0 ，a n dc o m p a r e di t s p e r f o r m a n c ew i t ht h o s eo ft h et r a d i t i o n a lG M R E S o nI n t e lC o r e2Q u a dC P UQ 9 4 0 0 2 6 6 G H za n d I n t e lC o r ei 7C P U9 2 0 2 6 7 G H z ，w h i c hs h o w e d4 0t i m e so fs p e e d - u pa n d2 0t i m e so fs p e e d - u po n a v e r a g er e s p e c t i v e l y K e yw o r d s ：C U D A ；G P G P U ；g e n e r a l i z e dm i n i m a lr e s i d u a lm e t h o d ；s p a r s em a t r i xv e c t o rm u l t i p l i c a t i o n 计算机图形学领域和一些实际工程应用中存在 很多偏微分方程的求解，如软体变形、流体仿真、几 何处理，这些方程在经过离散化后都转化成线性方 程组，从而将复杂问题求解变成一个可计算机求解 的问题这类求解通常采用迭代法进行数值计算，因 此此类方法的高效求解对这种复杂问题有着非常重 要的意义该线性方程组可统一表示为A x = b ；其中 A 为行X 恕大小的系数矩阵，z 为i t 元变量，b 为已知 收稿日期：2 0 1 0 0 9 2 5 ；修回日期：2 0 1 0 - 1 卜3 0 基金项目：国家自然科学基金( 6 0 9 7 3 0 6 6 。6 0 8 3 3 0 0 7 ) 中国科学院软件研究所计算机科学 国家重点实验宦开放基金( S Y S K F l 0 0 4 ) 柳有权( 1 9 7 6 一) ，男。博士，副教授C C F 会员，主要研究方向为计算机图形学、虚拟现实尹康学 ( 1 9 9 0 一) ，男，在校学生吴恩华( 1 9 4 7 一) ，男，博士，教授，博士生导师，C C F 高级会员，主要研究方向为计算机图形学、虚拟现实 万方数据 5 5 4计算机辅助设计与图形学学报第2 3 卷 量在目前诸多求解大规模线性方程组问题的迭代 法中，重开始广义极小残量法( g e n e r a l i z e dm i n i m a l r e s i d u a lm e t h o d ，G M R E S ) 1 1 是很受欢迎的算法 之一，它通过K r y l o v 子空间矢量的最小残量来迭代 求解，具有收敛速度快、稳定性好等优点 目前有很多研究人员侧重于改进G M R E S 算 法，以进一步提高该方法迭代的效率如全忠等 3 通 过构造多项式预处理因子来克服G M R E S 算法有 时收敛很慢或停滞的缺陷，H a b u 等 4 1 通过调整重 新开始来加快G M R E S 的收敛速度 但除了从算法层面来改进整体收敛速度，单步 计算的效率提升同样重要，这一方面依靠硬件的计 算能力的提升，另外新的计算架构通过对算法进行 并行化处理也可提升算法计算效率传统的高性能 计算依赖大型机和计算集群，然而这样的计算系统 都很昂贵G P U 的发展为高性能计算提供了另外一 种思路 5 弗 ，它采用众核架构，即芯片上集成了多个 并行处理单元随着可编程性的出现，G P U 从单纯的 图形流水线渲染转向通用计算上的应用( G P G P U ) ， 现在在高性能计算机领域也开始崭露头角G P U 单 位计算成本的下降引起了很多研究机构和企业的 广泛重视，如国产天河系统由于采用C P U + G P U 以较低的代价获得非常高的性能，在最近的全球高 性能计算机排行榜上名列第一N V I D I A L 7 0 推出的 C U D A ( c o m p u t eu n i f i e dd e v i c ea r c h i t e c t u r e ) 架构由 于编程方式的革新更是推动了G P U 芯片在高性能 计算上的应用 目前已有一些G M R E S 利用C U D A 进行加速 的工作，如W a n g 等在G e F o r c eG T X 2 8 0 图形卡上 获得2 0 倍加速 8 ，V e l a m p a r a m b i l 等类似的工作在 G e F o r c e8 8 0 0 上获得1 3 倍的加速L 9 ，G h a e m i a n 等 基于N V I D I AT e s l aC 8 7 0G P U 只获得6 0 的加 速 文献 1 1 对作为各种迭代求解方法都要使用 的核心算法稀疏矩阵与矢量乘法运算做了详细 分析虽然上述T 作各自的硬件平台略有不同，但整 体加速效率仍有提升的空间 本文基于C U D A 将G M R E S 算法在G P U 上进 行重新设计，尤其是对稀疏矩阵与矢量乘法部分，通 过合并访问和共享内存的分配来优化负载；并充分 利用G P U 的众核处理能力，使得算法效率相对于 C P U 算法有大幅度提升另外，通过跟N V I D I A 提 供的稀疏矩阵与矢量相乘的G P U 算法 1 1 比较，本 文算法实现相对于N V I D I A 自身的代码也有好几 倍的效率提升，同时代码在h t t p ：i m l a b c h d e d u c n c u d a 上公开，供研究人员免费使用 1G M R E S 算法分析 关于G M R E S 算法的详细描述请参考文献 2 ， 本文为了完整起见，对其稍作介绍对于线性方程组 A x = b ，G M R E S 算法的m 阶K r y l o v 子空间为k = s p a n ( b ，A b ，A 2 b ，A 一1 西) ；G M R E S 通过求使残量 A x 。一b 最小的矢量X 。K 。来逼近A x b 的精确 解但是，矢量b ，A b ，A 2 b ，A 一1b 几乎是线性相 关的，因此通常采用A r n o l d i 迭代方法来找出正交 矢量n ，耽，作为m 阶K r y l o v 子空间的基故 矢量工。k 可写成x 。= ，。Y 。，其中Y 。R ”且y 。 是由 ，1 ， ，2 ，1 ，。组成的行m 矩阵 通过A r n o l d i 迭代过程也可产生一个( m + 1 ) m 阶的上H e s s e n b e r g 矩阵H 。满足A V 。= L + 1H 。因 为，。是正交的，因此有| IA x 。b0 = I IH 。Y 。一 印。I I ；其中自一( 1 ，0 ，0 ，0 ) 是R 秆1 的标准基的 第一个矢量，并且口= l IA x 。一bI l ，x 。是初始矢量 ( 通常是零矢量) 因此，求使得残量，。= A x 。一b 范 数最小的工。，就变成求，。= H 。Y 。一胎。范数最小的 问题，即为一个m 阶线性最小二乘问题，通常情况 下仇咒 通过前面的描述可知，G M R E S 算法