摘 要： 偏微分方程数值解法（包括有限差分法、有限元法）以及大量的数学物理方程数值解法最终都会演变成求解大型线性方程组。因此，探讨快速、稳定、精确的大型线性方程组解法一直是数值计算领域不断深入研究的课题且具有特别重要的意义。在迭代法中，共轭斜量法（又称共轭梯度法）被公认为最好的方法之一。但是，该方法最大缺点是仅适用于线性方程组系数矩阵为对称正定矩阵的情况，而且常规的CPU算法实现非常耗时。为此，通过将线性方程组系数矩阵作转换成对称矩阵后实施基于GPU-CUDA的快速共轭斜量法来解决一般性大型线性方程组的求解问题。试验结果表明：在求解效率方面，基于GPU-CUDA的共轭斜量法运行效率高，当线性方程组阶数超过3000时，其加速比将超过14；在解的精确性与求解过程的稳定性方面，与高斯列主元消去法相当。基于GPU-CUDA的快速共轭斜量法是求解一般性大型线性方程组快速而非常有效的方法。
关键词： GPU； CUDA； 大型线性方程组； 共轭斜量法； 算法； 并行计算
中图分类号：TP311.1 文献标志码：A 文章编号：1006-8228（2014）04-04-03
Abstract： The numerical solution for partial differential equations （including finite difference method， and finite element method） and a large number of equations of mathematical physics problems will eventually evolve into solving a large-scale linear equation system. Therefore， studying fast， stable and accurate solutions for large-scale linear equation systems has been a hot topic in the field of numerical calculation for years， which has special significance. Among iterative methods， conjugate gradient method is recognized as one of the best methods. However， this method is only applicable to linear equation systems in which coefficient matrix is symmetric and positive definite. Besides， in conventional CPU implementation， the method for solving a large-scale linear equation system is time-consuming. After the linear equations coefficient matrix A is converted into a symmetric matrix by， the fast conjugate gradient method based on GPU-CUDA is implemented to solve a general large-scale linear equation system. The results show that it is highly efficient. When the rank of the coefficient matrix is over 3000， the speedup will be over 14 times. Besides， it has the same accuracy and stability as Gaussian elimination method. The conjugate gradient method based on GPU-CUDA becomes a fast and effective method for solving large-scale general linear equation systems.
Key words： GPU； CUDA； large-scale linear equation system； conjugate gradient method； algorithm； parallel computation
0 引言
偏微分方程数值解法（包括有限差分法、有限元法）及大量的数学物理方程数值解法最终都将演变成求解大型线性方程组[4，7]。随着离散网格数量的增加，线性方程组的阶数也同步增长。目前，求解实际问题的线性方程组的阶数一般都超过1000阶，因此，探讨稳定、快速、精确的大型线性方程组解法具有特别重要的意义。在迭代法中，共轭斜量法（亦称共轭梯度法）被公认为最好的方法之一。文献[5]等所述基于CPU的常规共轭斜量法已非常成熟，但是，该方法最大缺点是仅适用于线性方程组系数矩阵为对称正定矩阵的方程组，而且非常耗时，难以满足求解1000阶以上线性方程组的要求。在应用GPU新技术求解大型线性方程组方面，文献[6]成功应用GPU技术实现了LU分解和Laplace算法。文献[1]和文献[3]分别利用片元着色器编程和预处理平方根共轭梯度法（Preconditioned Conjugate Gradients Squared，PCGS），实现了基于GPU的共轭梯度法，有效加速了线性方程组的求解，但是，实现方法非常复杂。本文使用计算统一设备架构CUDA（Compute Unified Device Architecture）技术，将方程组系数矩阵A作ATA转换成对称矩阵后，通过CPU与GPU的协同来实现基于GPU-CUDA的快速共轭斜量法，实现方法简单并获得了较大的加速比。
1 共轭斜量法原理及算法[5，8]
共轭斜量法与最速下降法同属于一种方法—变分法。首先将线性方程组Ah=b的解h=A-1b（向量）看作是自变量x（向量）的某一个二次函数的极小点（使函数f（x）达到极小的点），然后用迭代法求f（x）的极小点，借以求出线性方程组Ah=b的解。
共轭斜量法是在最速下降法基础上发展起来的一种快速、有效的线性方程组迭代解法，并且是有限步收敛的迭代法。如下是共轭斜量法的推算过程：
根据上述共轭斜量法的算法理论，容易设计出基于CPU的串行算法程序。由于每次迭代都需要进行大量的矩阵与向量乘法、向量与向量的点积运算，所以，当线性方程组的阶数很大时，基于CPU的共轭斜量法非常耗时。
2 共轭斜量法的CUDA实现[1-2]
通过分析算法我们发现，迭代过程中大量的矩阵与向量乘法、向量与向量的点积运算都不需要使用计算完成的结果，适合于设计基于GPU-CUDA的并行算法，从而减少计算时间、提高算法的效率。根据共轭斜量法的推算公式以及CUDA程序设计思想，即CPU进行数据准备、初始化GPU设备和执行串行代码，GPU进行并行计算并向CPU返回计算结果。设计如下算法完成共轭斜量法基于GPU-CUDA的实现。
⑴ 主机端程序完成数据的初始化。
⑵ 申请三块与系数矩阵A相等容量的设备内存空间，分别用于存储系数矩阵A、系数矩阵的转置AT及变换结果ATA；申请六块与方程组右端向量b相等容量的设备内存空间，用于存储右端向量b、方程组的解向量x和每次迭代需要的中间向量。
⑶ 将主机端系数矩阵A及右端向量b拷贝回至设备端。
⑷ 设计内核并行程序，计算ATA和ATb。
⑸ 调用内核并行程序解线性方程组ATAx=ATb。
⑹ 将方程组的解向量x拷贝回主机。
基于GPU-CUDA的共轭斜量法的关键是并行处理网格的构造以及内核函数的设计，并行处理网格的构造及主机内核函数调用的主要源代码如下[6，10]。
其中，DIM为线性方程组的阶数。因篇幅所限，部分求解中间向量的内核函数源代码不再列出。算法实现过程中，为了使每次迭代过程都不进行主机与设备端的数据交换，提高总体并行处理性能，采取了如下三个并行优化措施：
⑴ 设备端一次性分配足够内存空间用于存储所有矩阵、右端向量以及中间向量；
⑵ 主机端按需要调用内核函数完成并行处理；
⑶ 并行处理过程中，主机与设备端不进行任何数据交换，一直等待所有并行处理结束后，才将线性方程组的解x从设备端拷贝回主机。
3 共轭斜量法的性能测试[9-10]
为了充分验证算法的正确性、稳定性，以及进行算法的性能测试，通过如下源代码随机生成任意阶线性方程组。
在基于Linux系统的Acer ASPIRE 4736G手提电脑（含512M内存、8个流处理器的Geforce G105M图形卡）上进行了大量共轭斜量法、高斯列主元消去法试验。试验结果表明：并行算法与串行算法求解线性方程组的结果一致，当线性方程组的阶数小于600时，基于GPU-CUDA的并行算法较CPU算法耗时；当阶数大于700后，并行算法处理效率明显高于串行算法，并具有较高的加速比。试验结果如表1所示。由于Geforce G105M仅是NVIDIA公司的低端GPU显卡，所以加速比不高，如果换成具有更多流处理器的GPU显卡，将会获得非常大的加速比。
通过与高斯列主元消去法的试算结果对比，得到如下重要结果。
⑴ 由于高斯列主元消去法中回代求解过程无法实现并行处理，基于GPU-CUDA的高斯列主元消去法比基于CPU的算法计算速度慢且稳定性差。因此，基于CPU的高斯列主元消去法仍然是目前解线性方程组最好的方法之一。
⑵ 对于共轭斜量法，当线性方程组阶数大于2000时，基于CPU的算法已相当耗时，基于GPU-CUDA的算法仍然具有6以上的加速比。线性方程组阶数越大，加速比越大，阶数超过3000阶时，加速比达14以上。
4 结束语
在分析共轭斜量法求解大型线性方程组原理及算法的基础上，设计并实现了基于GPU-CUDA的共轭斜量法的快速算法。通过大量的随机大型线性方程组的求解试验，得到以下两个具有重要意义的结论。
⑴ 基于GPU-CUDA的共轭斜量法不仅运行效率高，当线性方程组阶数超过700阶时，算法的加速比超过6，线性方程组阶数超过3000阶时，加速比达14以上。与此同时，基于GPU-CUDA的共轭斜量法在求解的精确性方面与列主元消去法相当。当线性方程组阶数大于3500后，尽管加速比变得更大，但是，其稳定性较CPU算法差，主要原因是基于GPU算法的误差积累所致。
⑵ 通过将一般性线性方程组系数矩阵A作ATA转换，右端项b作ATb转换后，即可应用基于GPU-CUDA的共轭斜量法求解线性方程组，从而拓展了共轭斜量法的应用范围，为求解一般性大型线性方程组提供了快速、有效的方法。
参考文献：
[1] NVIDIA. CUDA Programming Guide[OL].2012-12-22. http：///object/cuda_home.html.
[2] 夏健明，魏德敏著.共轭梯度法的GPU实现[J].微计算机工程，2009.17：274-276
[3] 李熙铭，欧阳丹彤，白洪涛著.基于GPU的混合精度平方根共轭梯度算法[J].仪器仪表学报，2012.1：98-104
[4] 李荣华，刘播著.微分方程数值解法（第四版）[M].高等教育出版社，2010.
[5] 蒋长锦编著.科学计算与C程序集[M].中国水利水电出版社，2010.
[6] 陈颖，林锦贤，吕暾著.LU分解和Laplace算法在GPU上的实现[J].计算机应用，2011.3：851-855
[7] 吕英华编著.计算电磁学的数值方法[M].人民教育出版社，2006.
[8] 武汉大学，山东大学著.计算方法[M].科学出版社，1983.
[9] 肖江，胡柯良，邓元勇著.基于CUDA的矩阵乘法和FFT性能测试[J].计算机工程，2009.10：7-10
[10] 刘丽，沈杰，李洪林著.基于GPU的矩阵求逆性能测试和分析[J].华东理工大学学报（自然科学版），2010.6：812-816