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

第２３卷第４期　２０１１年４月　计算机辅助设计与图形学学报　Ｊｏｕｒｎａｌ　ｏｆ　Ｃｏｍｐｕｔｅｒ～Ａｉｄｅｄ　Ｄｅｓｉｇｎ　８Ｌ　Ｃｏｍｐｕｔｅｒ　Ｇｒａｐｈｉｃｓ　Ｖｏ１．２３　ＮＯ．４　Ａｐｒ．２０１１　大规模稀疏线性方程组的ＧＭＲＥＳ—ＧＰＵ快速求解算法　柳有权　，尹康学”，吴恩华　。’　ｌ　（长安大学信息工程学院西安７１００６４）　北京１００１９０）　２　（中国科学院软件研究所计算机科学国家重点实验室’（澳门大学科技学院澳门）　（ｙｏｕｑｕａｎ＠ｃｈｄ．ｅｄｕ．ｃｎ）　摘要：重开始广义极小残量法（ＧＭＲＥｓ）是求解大规模线性方程组的常用算法之一，具有收敛速度快、稳定性好等　优点．文中基于ＣＵＤＡ将ＧＭＲＥＳ算法在ＧＰＵ上进行并行算法实现，尤其针对稀疏矩阵矢量乘法运算，通过合并访　问和共享内存策略相结合的手段使得算法效率大幅度提升．对于大规模数据集，在ＧｅＦｏｒｃｅ　ＧＴＸ　２６０上的运行结果　相对于Ｉｎｔｅｌ　Ｃｏｒｅ　２　Ｑｕａｄ　ＣＰＵ　Ｑ９４Ｏ０＠２．６６　ＧＨｚ得到了平均４Ｏ余倍的加速效果，相对于Ｉｎｔｅｌ　Ｃｏｒｅ　ｉ７　ＣＰＵ　９２０＠　２．６７　ＧＨｚ也可得到平均２Ｏ余倍的加速效果．　关键词：ＣＵＤＡ；ＧＰＧＰＵ；重开始广义极小残量法；稀疏矩阵矢量乘法　中图法分类号：ＴＰ３９１　Ｆａｓｔ　ＧＭＲＥＳ－ＧＰＵ　Ｓｏｌｖｅｒ　ｆｏｒ　Ｌａｒｇｅ　Ｓｃａｌｅ　Ｓｐａｒｓｅ　Ｌｉｎｅａｒ　Ｓｙｓｔｅｍｓ　Ｌｉｕ　Ｙｏｕｑｕａｎ　，　，Ｙｉｎ　Ｋａｎｇｘｕｅ”，ａｎｄ　Ｗｕ　Ｅｎｈｕａ　’。　（Ｓｃｈｏｏｌ　ｏｆ　Ｉｎｆｏｒｍａｔｉｏｎ　Ｅｎｇｉｎｅｅｒｉｎｇ，Ｃｈａｎｇ　ａｎ　Ｕｎｉｖｅｒｓｉｔｙ，ｘｉ＇ａｎ　７１００６４）　（Ｓｔａｔｅ　Ｋｅｙ　Ｌａｂｏｒａｔｏｒｙ　ｏｆ　Ｃｏｍｐｕｔｅｒ　Ｓｃｉｅｎｃｅ，Ｉｎｓｔｉｔｕｔｅ　ｏｆ　Ｓｏｆｔｗａｒｅ，Ｃｈｉｎｅｓｅ　Ａｃａｄｅｍｙ　ｏｆ　Ｓｃｉｅｎｃｅｓ，Ｂｅｉｊｉｎｇ　１００１９０）　。　（Ｆａｃｕｌｔｙ　ｆ　ｏＳｃｉｅｎｃｅ　ａｎｄ　Ｔｅｃｈｎｏｌｏｇｙ，Ｕｎｉｖｅｒｓｉｔｙ　ｏｆ　Ｍａｃａｕ，Ｍａｃａｏ）　Ａｂｓｔｒａｃｔ：Ａｓ　ａ　ｐｏｐｕｌａｒ　ｉｔｅｒａｔｉｖｅ　ｍｅｔｈｏｄ　ｔｏ　ｓｏｌｖｅ　ｌｉｎｅａｒ　ｅｑｕａｔｉｏｎｓ，ｒｅｓｔａｒｔｅｄ　ｇｅｎｅｒａｌｉｚｅｄ　ｍｉｎｉｍａｌ　ｒｅｓｉｄｕａｌ　ｍｅｔｈｏｄ（ＧＭＲＥＳ）ｈａｓ　ｔｈｅ　ａｄｖａｎｔａｇｅｓ　ｏｆ　ｆａｓｔ　ｃｏｎｖｅｒｇｅｎｃｅ　ａｎｄ　ｇｏｏｄ　ｓｔａｂｉｌｉｔｙ．Ｔｈｉｓ　ｐａｐｅｒ　ｉｍｐｌｅｍｅｎｔｓ　ａ　ｐａｒａｌｌｅｌ　ＧＭＲＥＳ　ｉｎ　ＧＰＵ　ｂａｓｅｄ　ｏｎ　ＣＵＤＡ．Ｐａｒｔｉｃｕｌａｒｌｙ，ｔｈｅ　ｓｐａｒｓｅ　ｍａｔｒｉｘ　ｖｅｃｔｏｒ　ｍｕｌｔｉｐｌｉｃａｔｉｏｎ　ｉｓ　ｏｐｔｉｍｉｚｅｄ　ｗｉｔｈ　ｃｏｈｅｒｅｎｃｅ　ｖｉｓｉｔｉｎｇ　ａｎｄ　ｓｈａｒｅｄ　ｍｅｍｏｒｙ，ｗｈｉｃｈ　ｓｉｇｎｉｆｉｃａｎｔｌｙ　ｉｍｐｒｏｖｅｓ　ｔｈｅ　ｐｅｒｆｏｒｍａｎｃｅ．Ｗｅ　ｔｅｓｔｅｄ　ｔｈｅ　ｐａｒａｌｌｅｌｅｄ　ＧＭＲＥＳ　ｏｎ　ａ　ＧＰＵ　ｏｆ　ＧｅＦｏｒｃｅ　ＧＴＸ２６０，ａｎｄ　ｃｏｍｐａｒｅｄ　ｉｔｓ　ｐｅｒｆｏｒｍａｎｃｅ　ｗｉｔｈ　ｔｈｏｓｅ　ｏｆ　ｔｈｅ　ｔｒａｄｉｔｉｏｎａ１　ＧＭＲＥＳ　ｏｎ　Ｉｎｔｅ１　Ｃｏｒｅ　２　Ｑｕａｄ　ＣＰＵ　Ｑ９４００＠２．６６ＧＨｚ　ａｎｄ　Ｉｎｔｅｌ　Ｃｏｒｅ　ｉ７　ＣＰＵ　９２０＠２．６７ＧＨｚ，ｗｈｉｃｈ　ｓｈｏｗｅｄ　４０　ｔｉｍｅｓ　ｏｆ　ｓｐｅｅｄ—ｕｐ　ａｎｄ　２０　ｔｉｍｅｓ　ｏｆ　ｓｐｅｅｄ—ｕｐ　ｏｎ　ａｖｅｒａｇｅ　ｒｅｓｐｅｃｔｉｖｅｌｙ．　Ｋｅｙ　ｗｏｒｄｓ：ＣＵＤＡ；ＧＰＧＰＵ；ｇｅｎｅｒａｌｉｚｅｄ　ｍｉｎｉｍａｌ　ｒｅｓｉｄｕａｌ　ｍｅｔｈｏｄ；ｓｐａｒｓｅ　ｍａｔｒｉｘ　ｖｅｃｔｏｒ　ｍｕｌｔｉｐｌｉｃａｔｉｏｎ　计算机图形学领域和一些实际工程应用中存在　很多偏微分方程的求解，如软体变形、流体仿真、几　的问题．这类求解通常采用迭代法进行数值计算，因　此此类方法的高效求解对这种复杂问题有着非常重　要的意义．该线性方程组可统一表示为Ａｘ＝ｂ；其中　Ａ为ｎ×　大小的系数矩阵，Ｘ为　元变量，ｂ为已知　何处理，这些方程在经过离散化后都转化成线性方　程组，从而将复杂问题求解变成一个可计算机求解　收稿日期：２０１０　０９　２５；修回日期：２０１０—１１－３０．基金项目：国家自然科学基金（６０９７３０６６，６０８３３００７）；中国科学院软件研究所计算机科学　国家重点实验室开放基金（ＳＹＳＫＦ１００４）．柳有权（１９７６一），男，博士，副教授，ＣＣＦ会员，主要研究方向为计算机图形学、虚拟现实．尹康学　（１９９ｏ一），男，在校学生．吴恩华（１９４７），男，博士，教授，博士生导师，ＣＣＦ高级会员，主要研究方向为计算机图形学、虚拟现实．　５５４　计算机辅助设计与图形学学报　第２３卷　量．在目前诸多求解大规模线性方程组问题的迭代　法中，重开始广义极小残量法（ｇｅｎｅｒａｌｉｚｅｄ　ｍｉｎｉｍａｌ　ｒｅｓｉｄｕａｌ　ｍｅｔｈｏｄ，ＧＭＲＥＳ）①＿】　是很受欢迎的算法　之一，它通过Ｋｒｙｌｏｖ子空间矢量的最小残量来迭代　求解，具有收敛速度快、稳定性好等优点．　目前有很多研究人员侧重于改进ＧＭＲＥＳ算　法，以进一步提高该方法迭代的效率．如全忠等口　通　过构造多项式预处理因子来克服ＧＭＲＥＳ算法有　时收敛很慢或停滞的缺陷，Ｈａｂｕ等＿４　通过调整重　新开始来加快ＧＭＲＥＳ的收敛速度．　但除了从算法层面来改进整体收敛速度，单步　计算的效率提升同样重要，这一方面依靠硬件的计　算能力的提升，另外新的计算架构通过对算法进行　并行化处理也可提升算法计算效率．传统的高性能　计算依赖大型机和计算集群，然而这样的计算系统　都很昂贵．ＧＰＵ的发展为高性能计算提供了另外一　种思路　］，它采用众核架构，即芯片上集成了多个　并行处理单元．随着可编程性的出现，ＧＰＵ从单纯的　图形流水线渲染转向通用计算上的应用（ＧＰＧＰＵ），　现在在高性能计算机领域也开始崭露头角．ＧＰＵ单　位计算成本的下降引起了很多研究机构和企业的　广泛重视，如国产天河系统由于采用ＣＰＵ＋ＧＰＵ　以较低的代价获得非常高的性能，在最近的全球高　性能计算机排行榜上名列第一．ＮＶＩＤＩＡ［７　推出的　ＣＵＤＡ（ｃｏｍｐｕｔｅ　ｕｎｉｆｉｅｄ　ｄｅｖｉｃｅ　ａｒｃｈｉｔｅｃｔｕｒｅ）架构由　于编程方式的革新更是推动了ＧＰＵ芯片在高性能　计算上的应用．　目前已有一些ＧＭＲＥＳ利用ＣＵＤＡ进行加速　的工作，如Ｗａｎｇ等在ＧｅＦｏｒｃｅ　ＧＴＸ２８０图形卡上　获得２０倍加速　］，Ｖｅｌａｍｐａｒａｍｂｉｌ等类似的工作在　ＧｅＦｏｒｃｅ　８８００上获得１３倍的加速ＬｇＪ，Ｇｈａｅｍｉａｎ等　基于ＮＶＩＤＩＡ　Ｔｅｓｌａ　Ｃ８７０　ＧＰＵ只获得６０　的加　速口　．文献Ｅｎ］对作为各种迭代求解方法都要使用　的核心算法——稀疏矩阵与矢量乘法运算做了详细　分析．虽然上述工作各自的硬件平台略有不同，但整　体加速效率仍有提升的空间．　本文基于ＣＵＤＡ将ＧＭＲＥＳ算法在ＧＰＵ上进　行重新设计，尤其是对稀疏矩阵与矢量乘法部分，通　过合并访问和共享内存的分配来优化负载；并充分　利用ＧＰＵ的众核处理能力，使得算法效率相对于　ＣＰＵ算法有大幅度提升．另外，通过跟ＮＶＩＤＩＡ提　供的稀疏矩阵与矢量相乘的ＧＰＵ算法［１妇比较，本　文算法实现相对于ＮＶＩＤＩＡ自身的代码也有好几　倍的效率提升，同时代码在ｈｔｔｐ：／／ｉｍｌａｂ．ｃｈｄ．ｅｄｕ．　ｃｎ／ｃｕｄａ／￣公开，供研究人员免费使用．　１　ＧＭＲＥＳ算法分析　关于ＧＭＲＥＳ算法的详细描述请参考文献［２］，　本文为了完整起见，对其稍作介绍．对于线性方程组　Ａｘ＝ｂ，ＧＭＲＥＳ算法的ｍ阶Ｋｒｙｌｏｖ子空间为Ｋ埘一　ｓｐａｎ（ｂ，Ａｂ，Ａ　ｂ，…，Ａ一　６）；ＧＭＲＥＳ通过求使残量　Ａｘ　一６最小的矢量Ｘ　∈　来逼近Ａｘ—ｂ的精确　解．但是，矢量ｂ，Ａｂ，Ａ　ｂ，…，Ａ一　６几乎是线性相　关的，因此通常采用Ａｒｎｏｌｄｉ迭代方法来找出正交　矢量＇，　，＇，２，…，ｖ　作为ｍ阶Ｋｒｙｌｏｖ子空间的基．故　矢量　∈Ｋ　可写成Ｘ　一Ｖ　Ｙ　，其中　∈　且ｙ　是由＇，１，’，２，…，＇，　组成的　×ｍ矩阵．　通过Ａｒｎｏｌｄｉ迭代过程也可产生一个（ｍ＋１）×ｍ　阶的上Ｈｅｓｓｅｎｂｅｒｇ矩阵Ｈ　满足ＡＶ　一’，珥＋１Ｈ　．因　为Ｖ　是正交的，因此有ｌｌ　Ａｘ　一ｂ　ｌＩ—ｌ　ｌＨ　Ｙ　一　ｌｌ；其中ｅ１一（１，０，０，…，ｏ）是Ｒ啪＋　的标准基的　第一个矢量，并且　一ｌｌ　Ａｘ。一ｂ　Ｉｌ，Ｘ。是初始矢量　（通常是零矢量）．因此，求使得残量ｒ　＝Ａｘ　一ｂ范　数最小的　，就变成求ｒ埘一Ｈ　Ｙ　一　范数最小的　问题，即为一个　阶线性最小二乘问题，通常情况　下ｍ《　．　通过前面的描述可知，ＧＭＲＥＳ算法在迭代的　每一步中要进行如下操作：　Ｓｔｅｐ１．做一步Ａｒｎｏｌｄｉ迭代计算（见算法１中的第②～　⑧步）．　Ｓｔｅｐ２．寻找使得ｌＩ　ｒ埘ｌｌ最小的Ｙ　（见算法１中的第⑨　步）．　Ｓｔｅｐ３．计算Ｘｍ—　０＋Ｖｍｙ　．　Ｓｔｅｐ４．如果残量不够小，增大Ｋｒｙｌｏｖ子空间维度并重　复以上步骤．　本文采用的是重开始ＧＭＲＥＳ算法，即将Ｓｔｅｐ４　改为：如果残量不够小，将　的初值置为　并重复　以上步骤．　算法１．ＧＭＲＥＳ算法　。　每次迭代的具体步骤　①计算ｒｏ＝６一Ａ　。，卢一ｌ　Ｉｒ０　ｌＩ和　一，。／卢；　②初始化（ｍ＋１）×ｍ阶的上Ｈｅｓｓｅｎｂｅｒｇ矩阵Ｈ　为　０；　③循环』一１，２，…，ｍ　④　计算　一Ａ　；　⑤　循环　一１，２，…，　①ｈｔｔｐ：／／ｅｎ．ｗｉｋｉｐｅｄｉａ．ｏｒｇ／ｗｉｋｉｌＧｅｎｅｒａｌｉｚｅｄ－ｍｉｎｉｍａｌ＿ｒｅｓｉｄｕａｌ＿ｍｅｔｈｏｄ　第４期　柳有权，等：大规模稀疏线性方程组的ＧＭＲＥＳ￣ＧＰＵ快速求解算法　⑥　计算ｈ　一　・ｖ　；计算　—ｗＪ—ｈｉｊｖ　⑦　计算ｈ　１．Ｊ＝ｌ１　ｌｌ，如果ｈ　一０，则将　赋给　ｍ，转到⑨；　⑧　计算　＋　—　／ｈ　⑨计算ｌ１日　Ｙ　一　ｌｌ的最小值，得到对应的　；　⑩计算　一Ｘ。＋Ｖ　．　设系数矩阵Ａ为ｎ阶可逆矩阵，每行非零元数　为ｋ，Ｋｒｙｌｏｖ子空间维度为ｍ（常数），则每次迭代中　各种计算所需次数及复杂度如表１所示．　表１　算法１中各计算步复杂度分析　很显然，无论从计算量还是从ＧＰＵ程序的优　化难度来说，表１中的前２项都是整个算法中最耗　时的部分．　设稀疏矩阵与矢量乘积、内积与范数的浮点运　算次数分别为厂　和厂　，则ｆ　一２ｎｋ（　＋１），ｆ　一　ｅｒ（ｍ＋３）＋２　．当忌》　一予＋１时，稀　疏矩阵与矢量乘积运算决定程序的整体性能；当　忌《　＋１时，矢量内积与范数运算决定程序的整体　性能．　考虑到矢量内积和范数运算在ＧＰＵ中运行比　稀疏矩阵和矢量乘积遇到的性能瓶颈要少很多，因　此可以认为ｍ值越大，加速效果越好；ｋ值越大，加　速效果越差．　尽管ｍ较大时可以加速收敛和提高精度，但同　时会使每次迭代的计算量增加并且使存储占用增　加；另外由于浮点精度有限，因此ｍ的值不宜取得　过大．至于ｍ的值取多少合适，要考虑矩阵本身的　特点、机器内存的、浮点精度的等，一般根　据矩阵条件数凭经验决定．另外，过小的ｍ可能导　致收敛停滞，文献Ｅ４］对此进行了一些讨论．　２　ＧＭＲＥＳ算法的ＣＵＤＡ优化实现　关于ＣＵＤＡ编程的详细说明可以参考ＮＶＩＤＩＡ　手册＿７　，本文不再赘述．根据本文前面的描述，ＧＭＲＥＳ　算法跟其他线性方程组迭代法求解类似．整个计算　被分解成稀疏矩阵与矢量乘积运算，矢量内积与范　数，稠密矩阵与矢量乘积，矢量加减、矢量乘标量，最　小二乘求解等几步．为了获得整体优化的效果，本文　对于不同的计算步骤采用不同的策略．　２．１　稀疏矩阵与矢量乘积运算　本文中矩阵存储使用压缩稀疏行存储（ｃｏｍｐｒｅｓｓｅｄ　ｓｐａｒｓｅ　ｒｏｗ，ＣＳＲ）格式，即将所有非零元素存入一　个数组ｄａｔａ，这些非零元素的列号存人另外一个数　组ｉｎｄｉｃｅｓ，而数组ｒｐｏｓ存放每一行第一个非零元素　在数组ｄａｔａ中的序号，ｒｐｏｓ中的最后一个元素存　放总的非零元素的个数．ＣＳＲ格式稀疏矩阵与矢量　乘积在ＧＰＵ中计算遇到的瓶颈以及一些算法实现　见文献［１１］．本文改进了翟艳堂等的稀疏矩阵与矢　量乘积运算方法①，对矩阵数据读取全部实现了合　并访问，取得了很好的加速效果．　乘法与加法分为两步的思想较为简单，其步骤　如下：　Ｓｔｅｐ１．Ｋｅｒｎｅｌ１．计算每个非零元与相应标量的乘积，保　存结果至全局存储器，用ｒａｒｅｓ表示．　Ｓｔｅｐ２．Ｋｅｒｎｅｌ２．读取Ｋｅｒｎｅｌｌ的结果，并将每行对应的　数据累加得到结果矢量．　Ｋｅｒｎｅｌ１较为简单，这里不再赘述，详细请参考　相关链接．对于Ｋｅｒｎｅｌ２，使用纹理存储并不能解决　非合并数据读取大量访问延时造成的性能瓶颈．　本文提出合并访问的Ｋｅｒｎｅｌ２算法（以ｂｌｏｃｋ大　小为１２８为例）．　算法２．稀疏矩阵与矢量乘积运算第２步核心　Ｓｔｅｐ１．ｔｉｄｂ＇￣－ｔｈｒｅａｄ　ｉｎｄｅｘ　ｉｎ　Ｂｌｏｃｋ　Ｓｔｅｐ２．ｔｉｄｇ￣ｔｈｒｅａｄ　ｉｎｄｅｘ　ｉｎ　Ｇｒｉｄ　Ｓｔｅｐ３．ｓｕｍ＂￣－－０．０　Ｓｔｅｐ４．ｄｅｆｉｎｅ　ｓｈａｒｅｄ　ａｒｒａｙ　ｒｐｏｓｓ［１２９］ａｎｄ　ｍｒｅｓ～Ｓ　ＥｌｚＣ　Ｓｔｅｐｓ．ｒｐｏｓ～ｓ［ｔｉｄ一６］・　—ｒ户０ｓ［ｔｉｄ—ｇ］　Ｓｔｅｐ６．ｉｆ　ｔｉｄｂ一０　ｔｈｅｎ　ｒｐｏｓ５［１２９］一ｒｐｏｓ［ｔｉｄ—ｇ＋　１２８］　Ｓｔｅｐ７．ｂａｓｅＡｄｄｒｅｓｓ￣－－ｒｐｏｓ５　Ｅｏ］　Ｓｔｅｐ８．ｒａｒｅｓ～５［ｔｉｄ一　一ｍｒ∞［ｂａｓｅＡｄｄｒｅｓｓ＋ｔｉｄ～６ｊ　Ｓｔｅｐ９．Ｐ一｛ｚ　Ｊ　ｒｆｌｏｓ＿ｓ［ｔｉｄ＿　≤ｘ＜ｒｐｏｓ一５［ｔｉｄ＿６＋１］）ｎ　｛３８ｌ　ｂａｓｅＡｄｄｒｅｓｓ￣ｘ￣ｂａｓｅＡｄｄｒｅｓｓ＋１２８）　Ｓｔｅｐ１０．ｉｆ　Ｐ一　ｇｏｔｏ　Ｓｔｅｐｌ３　①ｈｔｔｐ：／／　ｄａ　＿ｄｎ．　ｔ／ｃ０ｎｔ　。　ｐｒ０／　ｉｄｉａ＿　ｈ。　．　ｐ　。ｉｎｔｉｄ一５２　５５６　计算机辅助设计与图形学学报　第２３卷　Ｓｔｅｐｌ　１．ｆｏｒ　ａｌｌ　ｅｌｅｍｅｎｔｓ　ｉ　ｉｎ　Ｐ　约为５１个．假设第一个线程在ｂｌｏｃｋ中的编号为ｔｉｄ，　Ｓｔｅｐｌ２．　ｓ　ｍ＋一ｌｎｌ＇ｅｓｓ［　］　则ｔｉｄ　３２为其在ｗａｒｐ块中的编号．当ｔｉｄ　３２≤１２　Ｓｔｅｐｌ３．ｂａｓｅＡｄｄｒｅｓｓ￣－－ｂａｓｓＡｄｄｒｅｓｓ＋１２８　时，这５１个线程占用２个相邻ｗａｒｐ块；当ｔｉｄ　３２＞　Ｓｔｅｐｌ４．ｉｆ　ｂａｓｅＡｄｄｒｅｓｓ＜ｒｐｏｓ一５　Ｄｚ９３　ｇｏｔｏ　Ｓｔｅｐ８　１２时，这５１个线程占用３个相邻ｗａｒｐ块．因此，　Ｓｔｅｐｌ５．ｒｅｓｕｌｔ［ｔｉｄ—ｇ］一５Ｍｍ．　其中，ｒｐｏｓ—ｓ为每个ｂｌｏｃｋ分配的共享内存　Ｓｔｅｐｌ２中ＳＭ的利用率期望值为　１３＿（ｓｈａｒｅｄ　ｍｅｍｏｒｙ），用于存放该ｂｌｏｃｋ所处理的某些　３２　５１十，一１９　５１・一６４　１　３２　９６一Ｉ一●…．６４．　行第一个非零元素在数组ｄａｔａ中的序号；ｍｒｅｓ—Ｓ　类似地可以计算当ｋ一５时，ＳＭ的利用率约为　也是每个ｂｌｏｃｋ分配的共享内存，用于存放该ｂｌｏｃｋ　０．７７．　所处理的非零元与相应标量的乘积．通过共享内存，　另外，扩大ｍｒｅｓ—Ｓ的容量可以减少从Ｓｔｅｐ８～　可以大大加速数据访问，减少延迟．　Ｓｔｅｐｌ４的迭代次数，因此该算法效果很好．例如当　本文通过分析发现上述算法存在如下不足：假　ｍｒｅｓＳ容量为５１２，ｋ＜４时，从Ｓｔｅｐ８￣Ｓｔｅｐｌ４的迭　设稀疏矩阵平均每行非零元素为ｋ，则每个ｂｌｏｃｋ中　代在大多数情况下只需要进行一次．也就是说，所有　１　０ｏ　同时执行Ｓｔｅｐｌ２的线程约为　，对于是一１０，实际　需要的数据被一次性读入共享内存，然后所有线程　工作的只有１３个线程，远不足一个ｗａｒｐ（一个　同时对共享内存中的数据进行操作．　ｗａｒｐ由３２个线程组成）．也就是说从Ｓｔｅｐ８～　对Ｓｔｅｐ９求交集，可以根据｛　Ｉａ≤ｘ＜ｂ）ｎ｛ｚＩ　Ｓｔｅｐｌ４的每一次迭代中，当程序执行到Ｓｔｅｐｌ２时，　ｆ≤ｘ＜ｄ｝≠　甘一（６≤ｃ　Ｖ口≥　）很容易地完成．　１　０　Ｓｔｒｅａｍｉｎｇ　Ｍｕｈｉｐｒｏｃｅｓｓｏｒｓ（ＳＭ）的利用率只有　．　本文与文献［１１］中的稀疏矩阵和矢量乘积算法　进行了代码级性能比较，结果如表２所示，其中本文　更糟糕的是，Ｓｔｅｐｌ２是上述算法计算量最多的一步．　所用的矩阵数据全部来自ｈｔｔｐ：／／ｗｗｗ．ｃｉｓｅ．ｕｆ１．　解决的办法是扩大ｍｒｅｓ—ｓ的容量．例如ｍｒｅｓ—Ｓ　ｅｄｕ／ｒｅｓｅａｒｃｈ／ｓｐａｒｓｅ／ｍａｔｒｉｃｅｓ．　的长度为５１２，对于是一１０，同时执行Ｓｔｅｐ１２的线程数　表２稀疏矩阵和矢量乘积算法效率比较　由表２可以看出，本文实现的稀疏矩阵和矢量乘　况下将每个标量乘积累加起来．本文采用ｒｅｄｕｃｔｉｏｎ￣　积算法要比ＮＶＩＤＩＡ提供的算法效率高，而且对于一　的思想，即将２个输入矢量划分成若干对小矢量，每　些矩阵，ＮⅥＤＩＡ提供的算法存在失效的情况．　个ｂｌｏｃｋ负责计算一对小矢量的内积，这些小矢量　２．２矢量内积和范数运算　的内积被写到ｍａｐｐｅｄ　ｍｅｍｏｒｙ中，由ＣＰＵ负责将　算法中矢量范数采用欧几里德范数，由于它与矢　这些小矢量内积加起来得到最终结果．　量内积计算类似，因此本文只介绍矢量内积的计算．　内积计算的困难在于如何在保持一定并行度的情　第４期　柳有权，等：大规模稀疏线性方程组的ＧＭＲＥＳ－ＧＰＵ快速求解算法　为了对不同问题规模都能得到最佳的ｗａｒｐ块　装载量，本文并未固定小矢量的长度，而是在　ｂｌｏｃｋ个数不大于５１２的情况下动态地决定小矢量　设输入矢量长度为ｎ，ｂｌｏｃｋ大小为ｂ，ｂｌｏｃｋ　个数为ｇ，小矢量长度为ｍ，则有如图１所示的判　定树．　的长度，因此得到的小矢量不多于５１２对．　图１矢量内积计算判定树　从图１中可以看出，利用该判定树得到的线程　和数据划分在任何时候都不会对ｗａｒｐ块装载量产　生影响．当，２≥５１２×５１２时，尽管此时ｂｌｏｃｋ的数量　误差舍入．例如，有一个长度为１　０００　０００的数组，每　个元素均为０．００１，使用ＩＥＥＥ标准单精度顺序相加　得到的结果是９９１．１４１　５４１，和正确结果１　０００相差　被为５１２，但由于每个ＳＭ只能同时装载２个　大小为５１２的ｂｌｏｃｋ，所以这个是合理的，不会　对性能产生影响．　接近百分之一．这是由于浮点精度导致的舍入误差　造成的．使用ｒｅｄｕｃｔｉｏｎ方法的情况要好很多，如算　法３代码所示，数组元素是两两相加的．当数组元素　的数量级差别不大时，使用ｒｅｄｕｃｔｉｏｎ方法得到的结　果几乎是精确的；但当数组元素数量级相差较大时，　其结果是偏小的，如１．ＯＥＩＯ＋１．ＯＥ一１０—１．０Ｅｌ０．解　以ｂｌｏｃｋ大小为１２８为例（图１判定树根的右　节点），ｋｅｒｎｅｌ函数代码如下：　算法３．ＣＵＤＡ内积运算核心　ｇｌｏｂａｌｖｏｉｄ　ｉｎｎｅｒＰｒｏ１２８（ｆｌｏａｔ＊ｖ，ｆｌｏａｔ＊＇‘，，ｆｌｏａｔ＊　决加法误差的方法有Ｋａｈａｎ　ｓｕｍｍａｔｉｏｎ算法①，但　其会造成性能上的大幅度下降，因此本文并未采用．　２．３稠密矩阵矢量乘积、矢量加减、矢量乘标量和　ｔｌ￣ｔｒｅ￥）｛　ｉｎｔ　ｔｉｄ—ｔｈｒｅａｄｌｄｘ．ｚ；　ｉｎｔ　ｔｉｄｉｎｇｒｉｄ—ｂｌｏｃｋＩｄｘ．ｚ＊ｂｌｏｃｋＤｉｍ．ｚ＋　最小二乘求解　ｔｈｒｅａｄｌｄｘ．ｚ：　ｓｈａｒｅｄｆｌｏａｔ　ｒ　ｓ［１２８］；　稠密矩阵矢量乘积、矢量加减和矢量乘标量问　题较为简单，本文采用ＣＵＤＡ实现．ＣＵＤＡ　ＳＤＫ中　提供了类似的例子代码，这里不在赘述．求解最小二　ｒ　ｓ［ｔｉｄ］一ｖ［ｔｉｄｉｎｇｒｉｄ￣＊ｗ［ｔｉｄ—ｉｎ—ｇｒｉｄ］；　ｓｙｎｃｔｈｒｅａｄｓ（）；　ｉｆ（ｔｉｄ＜６４）ｒＳ［ｔｉｄ］＋一ｒ—Ｓ［ｔｉｄ＋６４］；　乘问题Ｙ—ａｒｇ　ｒａｉｎ　Ｉ　ｌＨｍｙ—ｔ３ｅ　ｌ　较为烦琐口　］，ｌ　但从表１中可以看出，最小二乘问题的浮点运算次　数大约为４ｍ。，计算量并不大，没有构成问题的瓶　ｓｙｎｃｔｈｒｅａｄｓ（）；　ｉｆ（ｔｉｄ＜３２）｛　ｒ　ｓ［ｔｉｄ］＋一ｒ＿５［￡　＋３２］；　ｉｆ（ｔｉｄ＜１６）ｒｓ［ｔｉｄ］＋一ｒ—ｓ［￡　＋１６３；　颈，因此该计算在ＣＰＵ端完成．　ｉｆ（ｔｉｄ＜８）ｒｓ［ｔｉｄ］＋一ｒ＿＿ｓ［￡　＋８］；　ｉｆ（ｔｉｄ＜４）ｒｓ［ｔｉｄ］＋一ｒ－ｓ［ｆ　＋４］；　３实验与结果讨论　本文算法在Ｖｉｓｕａｌ　Ｓｔｕｄｉｏ　２００５＋Ｗｉｎｄｏｗ　ｉｆ（ｔｉｄ＜２）ｒｓＥｔｉｄ］＋一ｒ—ｓ［￡　＋２］；　ｉｆ（ｔｉｄ￣１）ｍｒｅｓ［ｂｌｏｃｋＩｄｘ．ｚ］一ｒ一　［ｏ］＋ｒ—ｓｉｌｌ；　ＸＰ／Ｖｉｓｔａ软件环境下，采用Ｃ＋＋语言结合ＣＵＤＡ　）　３．１编程实现．ＧＰＵ为ＮＶＩＤＩＡ　ＧｅＦｏｒｃｅ　ＧＴＸ２６０，　分别在２台装有不同ＣＰＵ的机器上运行，机器１为　其中　和Ｗ分别为输入的２个矢量，而ｍｒｅｓ则　为生成的小矢量．　矢量内积计算遇到的另一个问题是数组求和的　５５８　计算机辅助设计与图形学学报　第２３卷　Ｉｎｔｅｌ　Ｃｏｒｅ　２　Ｑｕａｄ　ＣＰＵ　Ｑ９４００＠２．６６　ＧＨｚ，得到了　ＣＰＵ　９２０＠２．６７　ＧＨｚ，得到平均２０余倍的加速效　平均４０余倍的加速效果，机器２为Ｉｎｔｅｌ　Ｃｏｒｅ　ｉ７　果．它们均迭代１Ｏ次，具体实验结果如表３所示．　表３机器１。２上的性能测试　０．２３４　０．５６３　２１．４０７　０．６７２　６１．７８２　３．０６２　３２　２０８．Ｏ３１　４．７８２　４３．５　由于机器２的ＣＰＵ配置比机器１高１倍，导致　从表４可以看出，曲线收敛速度跟矩阵关系比　机器２上的ＧＭＲＥＳ算法的ＣＰＵ版本运行的效率　较紧密，有些收敛得快，有些收敛得慢，说明收敛效　比机器１的提升了２倍，但同样的ＧＰＵ版本代码效　果首先取决于方程组系数矩阵的实际情况．对于同　率却略有下降，这可能是数据总线传输速率差异导　一个矩阵来说，Ｋｒｙｌｏｖ子空间维度的大小则对于收　致的．　敛有着至关重要的影响．另外可以看出，ｒａｊａｔ２６和　为了确保ＧＰＵ与ＣＰＵ运算结果一致，本文进　ｃ一２６　２个矩阵的求解并没有收敛到近似解，这是矩　行了收敛曲线和误差分布分析，总共迭代１００次，取　阵的条件数较大对误差较为敏感造成的，从收敛曲线　了５个矩阵样本；其中测试矩阵得到的收敛曲线和　也可看到这２个矩阵的计算收敛情况出现异常．　误差分步直方图如表４所示．　第４期　柳有权，等：大规模稀疏线性方程组的ＧＭＲＥＳ－ＧＰＵ快速求解算法　５５９　表４收敛曲线与误差分步直方图　收敛曲线　（横坐标为迭代次数，０　纵坐标为ｌｇ＿＿ｃＰ　　子空间维度　）　误差分步直方图　Ｏ　（横坐标为标量误差Ｉ　　Ｉ。．　吨如　ｌ，纵坐标为所占比例）　１．Ｏ　Ｏ．８　ｌ　ｌ　ｌ　０．６　０．４　ｌ　ｌ　０．２　｝　＼　１　１１　２１　３１　４１　５１　６１　７１　８１　９１　Ｌ　Ｉ　Ｏ　≥ｌ＿０　｜＿　Ｊ　ｌ　ｌ—Ｌ＿ｌＩ＿－ＩＪ—■　１．０Ｅ一２　１．０１３－４　１．０Ｅ一６≤１．０Ｅ一７　Ｉ－　■　Ｉ　　　．１．Ｏ　Ｏ．８　１　３　５　７　ｋ—～　、～　——～—０．６　１６　—　——　０．４　．　ｌ　９　｜　、　＼　＿　Ｏ．２　ｌ　ｌ　≥１．０　Ｉ　＿　１　．　１　１．０Ｅ一２　１．０Ｅ一４　１．０Ｅ一６《１．０Ｅ一７　１１　Ｏ　１１　２１　３１　４ｌ　５１　６１　７１　８１　９１　Ｌ　Ｊ　ｌ　ｌ　－　Ｌ一■．１　１　ｌ　ｌ　Ｏ．８　Ｏ．６　’●Ｌ　、　１．Ｏ　～ｌ　０．４　３２　—　＾＾　Ｎ～．　ｌ　　０．２　ｌＯ　ｌ　Ｉ　ｌ　Ｌ　．—，－＿、，＿　＾　一　ｌ一　≥１．０　．　Ｉ　ｌ　Ｌ　＿　一　１．０Ｅ一２　１．０１３—４　１．０Ｅ一６《１．０Ｅ一７　１　１１　２１　３１　４１　５１　６１　７１　８１　９１　注：——ｃ－２６；——ｂｃｓｓｔｋｌ５；－Ｃ－２６：　＿ｂｃｓｓｔｋｌ５：　——ｃｒｙｓｔｋ０２；——ｒａｊａｔ２６；——ＡＳＩＣ＿６８０ｋｓ；　’ｃｒｙｓｔｋ０２；　－ｒａｊａｔ２６；　－ＡＳＩＣ＿６８０ｋｓ　我们拟将本文中的求解算法应用到流体仿真和变形　４结论和未来工作　本文通过对ＧＭＲＥＳ算法进行算法分析，并利　用ＣＵＤＡ实现线性方程组的快速求解．尤其是对于　稀疏矩阵与矢量乘积运算部分，通过合并访问和共　享内存策略提高数据负载的利用率，使得对于大规　模线性方程组的求解效率得到提升，对于当前主流　的ＰＣ机来说，ＧＰＵ版本要比ＣＰＵ版本快２０～４０　倍，对于某些矩阵效率提升更多．由于ＧＭＲＥＳ迭　计算的计算机动画研究中去．　致谢感谢ＮＶＩＤＩＡ公司提供ＧＴＸ２６０显卡，　感谢佛罗里达大学提供大规模稀疏矩阵测试数据！　参考文献（Ｒｅｆｅｒｅｎｃｅｓ）：　Ｓａａｄ　Ｙ，［１］　Ｓｃｈｕｈｚ　Ｍ　Ｈ．ＧＭＲＥＳ：ａ　ｇｅｎｅｒａｌｉｚｅｄ　ｍｉｎｉｍａｌ　ｒｅｓｉｄｕａｌ　ａｌｇｏｒｉｔｈｍ　ｆｏｒ　ｓｏｌｖｉｎｇ　ｎｏｎｓｙｍｍｅｔｒｉｃ　ｌｉｎｅａｒ　ｓｙｓｔｅｍｓ　［Ｊ］．ＳＩＡＭ　Ｊｏｕｒｎａｌ　ｏｎ　Ｓｃｉｅｎｔｉｆｉｃ　ａｎｄ　Ｓｔａｔｉｓｔｉｃａｌ　Ｃｏｍｐｕｔｉｎｇ，　１９８６，７（３）：８５６－８６９　代算法适用于一般的线性方程组的求解，所以只要　求系数矩阵可逆就有解．虽然条件数对迭代收敛结　ｅｒａｔｉｖｅ　ｍｅｔｈｏｄｓ　ｆｏｒ　ｓｐａｒｓｅ　ｌｉｎｅａｒ　ｓｙｓｔｅｍｓ［Ｍ］．２ｎｄ　Ｅ２］　Ｓａａｄ　Ｙ．Ｉｔｅｄ．Ｐｈｉｌａｄｅｌｐｈｉａ：ＳＩＡＭ，２００３　果会有影响，相信在工程应用领域会大有用武之地．　５６Ｏ　计算机辅助设计与图形学学报　第２３卷　１－３３　Ｑｕａｎ　Ｚｈｏｎｇ，Ｘｉａｎｇ　Ｓｈｕｈｕａｎｇ．Ａ　ＧＭＲＥＳ　ｂａｓｅｄ　ｐｏｌｙｎｏｍｉａｌ　ｐｒｅｃ０ｎｄｉｔｉ０ｎｉｎｇ　ａｌｇｏｒｉｔｈｍ［Ｊ］．Ｍａｔｈｅｍａｔｉｃａ　Ｎｕｍｅｒｉｃａ　Ｓｉｎｉｃａ，　２００６，２８（４）：３６５—３７６（ｉｎ　Ｃｈｉｎｅｓｅ）　（全忠，向淑晃．基于ＧＭＲＥＳ的多项式预处理广义极小残　差法［Ｊ３．计算数学，２００６，２８（４）：３６５—３７６）　［４］Ｈａｂｕ　Ｍ，Ｎｏｄｅｒａ　Ｔ．ＧＭＲＥＳ（ｍ）ａｌｇｏｒｉｔｈｍ　ｗｉｔｈ　ｃｈａｎｇｉｎｇ　ｔｈｅ　ｒｅｓｔａｒｔ　ｃｙｃｌｅ　ａｄａｐｔｉｖｅｌｙ　Ｅｃ］／／Ｐｒｏｃｅｅｄｉｎｇｓ　ｏｆ　Ａｌｇｏｒｉｔｍｙ　Ｃｏｎｆｅｒｅｎｃｅ　ｏｎ　Ｓｃｉｅｎｔｉｆｉｃ　Ｃｏｍｐｕｔｉｎｇ．Ｈｅｉｄｅｌｂｅｒｇ：Ｓｐｒｉｎｇｅｒ，　２０００：２５４－２６３　［５３　Ｗｕ　Ｅｎｈｕａ，Ｌｉｕ　Ｙｏｕｑｕａｎ．Ｇｅｎｅｒａｌ　ｐｕｒｐｏｓｅ　ｃｏｍｐｕｔａｔｉｏｎ　ｏｎ　ＧＰＵ［Ｊ］．Ｊｏｕｒｎａｌ　ｏｆ　Ｃｏｍｐｕｔｅｒ—Ａｉｄｅｄ　Ｄｅｓｉｇｎ＆Ｃｏｍｐｕｔｅｒ　Ｇｒａｐｈｉｃｓ，２００４，１６（５）：６０１—６１２（ｉｎ　Ｃｈｉｎｅｓｅ）　（吴恩华，柳有权．基于图形处理器（ＧＰＵ）的通用计算［Ｊ］．　计算机辅助设计与图形学学报，２００４，１６（５）：６０１—６１２）　［６］Ｗｕ　Ｅ　Ｈ，Ｌｉｕ　Ｙ　Ｑ．Ｅｍｅｒｇｉｎｇ　ｔｅｃｈｎｏｌｏｇｙ　ａｂｏｕｔ　ＧＰＧＰＵ　［Ｃ］／／Ｐｒｏｃｅｅｄｉｎｇｓ　ｏｆ　ＩＥＥＥ　Ａｓｉａ　Ｐａｃｉｆｉｃ　Ｃｏｎｆｅｒｅｎｃｅ　ｏｎ　Ｃｉｒｃｕｉｔｓ　ａｎｄ　Ｓｙｓｔｅｍｓ．Ｌｏｓ　Ａｌａｍｉｔｏｓ：ＩＥＥＥ　Ｃｏｍｐｕｔｅｒ　Ｓｏｃｉｅｔｙ　Ｐｒｅｓｓ，　２００８：６１８—６２２　［７］ＮＶＩＤＩＡ　ＣＵＤＡ　Ｃ　ｐｒｏｇｒａｍｍｉｎｇ　ｇｕｉｄｅ．Ｖｅｒｓｉｏｎ　３．１［Ｍ］．　Ｓａｎ　Ｊｏｓｅ：ＮＶＩＤＩＡ，２０１０　［８］Ｗａｎｇ　Ｍ　Ｌ，Ｋｌｉｅ　Ｈ，Ｐａｒａｓｈａｒ　Ｍ，ｅｔ　ａ１．Ｓｏｌｖｉｎｇ　ｓｐａｒｓｅ　ｌｉｎｅａｒ　ｓｙｓｔｅｍｓ　ｏｎ　ＮＶＩＤＩＡ　ｔｅｓｌａ　ＧＰＵｓ［Ｍ］／／Ｌｅｃｔｕｒｅ　Ｎｏｔｅｓ　ｉｎ　Ｃｏｍｐｕｔｅｒ　Ｓｃｉｅｎｃｅ．Ｈｅｉｄｅｌｂｅｒｇ：Ｓｐｒｉｎｇｅｒ，２００９，５５４４：８６４—　８７３　［９］Ｖｅｌａｍｐａｒａｍｂｉｌ　Ｓ，ＭａｃＫｉｎｎｏｎ－Ｃｏｒｍｉｅｒ　Ｓ，Ｐｅｒｒｙ　Ｊ，ｅｔ　ａ１．　ＧＰＵ　ａｃｃｅｌｅｒａｔｅｄ　Ｋｒｙｌｏｖ　ｓｕｂｓｐａｃｅ　ｍｅｔｈｏｄｓ　ｆｏｒ　ｃｏｍｐｕｔａｔｉｏｎａｌ　ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃｓ［ｃ］／／Ｐｒｏｃｅｅｄｉｎｇｓ　ｏｆ　ｔｈｅ　３８ｔｈ　Ｅｕｒｏｐｅａｎ　Ｍｉｃｒｏｗａｖｅ　Ｃｏｎｆｅｒｅｎｃｅ．Ｌｏｓ　Ａｌａｍｉｔｏｓ：ＩＥＥＥ　Ｃｏｍｐｕｔｅｒ　Ｓｏｃｉｅｔｙ　Ｐｒｅｓｓ，２００８：１３１２—１３１４　［１０］Ｇｈａｅｍｉａｎ　Ｎ，Ａｂｄｏｌｌａｈｚａｄｅｈ　Ａ，Ｈｅｉｎｅｍａｎｎ　Ｚ，ｅｔ　ａ１．　Ａｃｃｅｌｅｒａｔｉｎｇ　ｔｈｅ　ＧＭＲＥＳ　ｉｔｅｒａｔｉｖｅ　ｌｉｎｅａｒ　ｓｏｌｖｅｒ　ｏｆ　ａｎ　ｏｉｌ　ｒｅｓｅｒｖｏｉｒ　ｓｉｍｕｌａｔｏｒ　ｕｓｉｎｇ　ｔｈｅ　ｍｕｌｔｉ—ｐｒｏｃｅｓｓｉｎｇ　ｐｏｗｅｒ　ｏｆ　ｃｏｍｐｕｔｅ　ｕｎｉｆｉｅｄ　ｄｅｖｉｃｅ　ａｒｃｈｉｔｅｃｔｕｒｅ　ｏｆ　ｇｒａｐｈｉｃｓ　ｃａｒｄｓ［ｃ］／／　Ｐｒｏｃｅｅｄｉｎｇｓ　ｏｆ　ｔｈｅ　９ｔｈ　Ｉｎｔｅｒｎａｔｉｏｎａｌ　Ｗｏｒｋｓｈｏｐ　ｏｎ　Ｓｔａｔｅ—ｏｆ—ｔｈｅ—Ａｒｔ　ｉｎ　Ｓｃｉｅｎｔｉｆｉｃ　ａｎｄ　Ｐａｒａｌｌｅｌ　Ｃｏｍｐｕｔｉｎｇ．　Ｈｅｉｄｅｌｂｅｒｇ：Ｓｐｒｉｎｇｅｒ，２００８：１５６—１５９　［１１］Ｂｅｌｌ　Ｎ，　Ｇａｒｌａｎｄ　Ｍ．　Ｅｆｆｉｃｉｅｎｔ　ｓｐａｒｓｅ　ｍａｔｒｉｘ—ｖｅｃｔｏｒ　ｍｕｌｔｉｐｌｉｃａｔｉｏｎ　ｏｎ　ＣＵＤＡ　ＪＲ］．　Ｓａｎ　Ｊｏｓｅ：　ＮＶＩＤＩＡ，　ＮＶＲ一２００８—００４。２００８　［１２］ＮＶＩＤＩＡ　ＣＵＤＡ　Ｃ　ｂｅｓｔ　ｐｒａｃｔｉｃｅｓ　ｇｕｉｄｅ．Ｖｅｒｓｉｏｎ　３．１［Ｍ］　Ｓａｎ　Ｊｏｓｅ：ＮＶＩＤＩＡ，２０１０