, 王文琛, 汪宇航, 王瑜华北电力大学 电子与通信工程系, 保定 071001
收稿日期：2019-09-08
基金项目：河北省自然科学基金资助项目（F2017502043）；计算机体系结构国家重点实验室开放课题（CARCH201802）；中央高校基本科研业务费专项资金项目（2017MS114）
作者简介：刘珊瑕(1994—), 女, 硕士研究生
通讯作者：靳松, 副教授, E-mail:jinsong@ncepu.edu.cn
摘要：稀疏线性方程组的求解是许多大规模科学计算任务的核心环节。目前，并行算法的发展为稀疏线性方程组的求解提供了新的思路和强有力的工具。然而，现有的并行算法存在一些缺陷，如最优子矩阵的划分难以获得、并行任务间的同步开销较大等。针对上述问题，该文提出一种基于变量相关性分解方法的稀疏线性方程组并行求解算法。该算法首先对系数矩阵进行不完全LU分解，得到上三角和下三角方程组，然后在这2个方程组求解过程中利用y与x的关系分解变量的相关性，同时并行计算变量的独立部分值，最后将所有的独立部分值相加得到变量的最终值。由于算法中变量的求解无需等待其所有前继变量计算完成即可进行部分值计算，因此有效减少了算法的执行时间，进而提高了算法的求解速度及并行度。实验结果表明：与调用cusparse库函数实现的并行求解方法相比，该文提出的算法能将稀疏线性方程组的求解速度提升了50%以上。
关键词：图形处理器稀疏线性方程组并行计算变量部分值线程
Parallel algorithm for solving sparse linear equations based on variable correlation decomposition
LIU Shanxia, JIN Song
, WANG Wenchen, WANG Yuhang, WANG YuDepartment of Electronics and Communication Engineering, North China Electric Power University, Baoding 071001, China
Abstract: Sparse linear equations are the core of many large scientific computing tasks. Although parallel algorithms are providing powerful, new tools for solving sparse linear equations, existing parallel algorithms have some drawbacks such as the optimal sub-matrix being difficult to obtain and the algorithms requiring a large overhead to synchronize parallel tasks. This paper presents a parallel algorithm for sparse linear equations based on the variable correlation decomposition method. Firstly, the algorithm performs an incomplete LU decomposition on the coefficient matrix to obtain two equations for the upper and lower triangles. Then, the correlation between y and x is used to solve for the independent part of the variable in parallel from the upper and lower triangles. All the individual partial values are then added to get the final value of the variable. Since the solution does not need to wait for all the previous variables to be calculated, the partial value calculation can be performed in parallel as the calculation proceeds, which significantly reduces the algorithm execution time and improves the algorithm solution speed and parallelism. Tests show that this sparse linear equation solver is more than 50% faster than the parallel solution method in the cusparse library.
Key words: graphics processing unit (GPU)sparse linear equationsparallel computationvariable partial valuethread
大规模稀疏线性方程组的求解是许多科学和工程计算中最重要的也是最基本的问题之一[1]。在一些以预处理迭代[2]算法为基础的科学计算任务中，方程组求解时间往往在求解整个问题所需的时间中占据很大的比重，有时甚至高达80%以上，成为提高计算速度的瓶颈[3]。由此可见，提升稀疏线性方程组的求解速度对优化整个科学计算任务的性能至关重要。目前，如何快速高效地求解稀疏线性方程组成为研究的热点[4]。
在求解稀疏线性方程组的硬件计算平台方面，以图形处理器(graphics processing unit，GPU)[5]为代表的众核处理器[6]在通用计算领域[7]展示了强大的优势，其丰富的并行计算资源为求解大规模稀疏线性方程组提供了强有力的支撑。然而，已有的计算方法，无论是来自学术界还是工业界，均受限于方程组求解的传统观点，即某一个变量的求解必须等到其所有前继变量求解完成之后才能开始。这种方法一定程度上限制了求解稀疏线性方程组的任务并行度，导致无法充分利用众核处理器丰富的并行硬件资源。为了解决变量求解时依赖关系而导致的顺序问题，现有的研究大多数集中在添加预处理阶段，将变量的求解划分为若干个集合(称为水平集(level-set)[8-9]或颜色集(color-set)[10-11])。即使必须按顺序执行集合，也可以并行计算任何单个集合中的变量，进而可以有效地利用并行硬件资源。这类方法通常比CPU[12-13]和GPU[14-15]上按照变量依赖关系的顺序求解要好得多。然而，基于集合的方法存在2个性能问题：首先，获得最优子矩阵的划分通常需要花费太多的时间，这有可能会抵消甚至消除并行计算的优势；其次，水平集之间的同步[16-17]开销降低了算法的并行化效率。
针对已有方法存在的问题，本文提出一种基于变量相关性分解方法的稀疏线性方程组并行求解算法。该算法通过分解变量求解的依赖关系，把求解稀疏线性方程组的过程划分成多个独立的计算任务，多个计算任务之间并行执行。该算法由于变量计算时无需等待其所有前继变量计算完成，而且不需要分析矩阵的复杂结构，完全消除了预处理阶段，有效减少了算法的执行时间，进而提高了算法的计算速度和并行度。本文在GPU实验平台实现了提出的并行算法，基于Florida稀疏矩阵集[18]构造了稀疏线性方程组，并应用提出的算法进行求解。实验结果表明：与调用cusparse[19]库函数实现的并行求解方法相比，本文提出的算法能将计算速度提升50%以上。
1 研究背景和现状在科学计算和工程技术中，许多问题的解决最终都转化为大规模稀疏线性系统的求解，如计算流体力学、计算电磁学、最优化问题、线性弹性力学、数值天气预报及核爆数值模拟等[20]。基本求解方法分为2种:迭代法[21-22]和直接法[23-24]。
迭代法一般通过构造迭代公式逐次逼近线性方程组的理论解得到近似解。迭代法求解方程组时通常与预处理技术结合。如针对多GPU平台上的Possion问题，Ament等[25]提出了具有不完全Possion预处理的共轭梯度(preconditioned conjugate gradient，PCG)方法，研究了采用块不完全LU分解预条件广义最小残量(generalized minimum residual，GMRES)[26-27]方法求解面向GPU平台的稀疏线性系统；Dziekonski等[28]给出了在有限元分析中应用Jacobi和Gauss-Seidel法对复共轭梯度进行预条件处理的方法。然而，迭代次数和收敛速度以及近似解与理论解的偏差等问题会影响算法的计算精度和速度；Gaikwad等[29]提出了一种基于Krylov子空间迭代求解器的实现方法，用于求解由该方法产生的多个小型系统，但是通信的开销会影响求解的性能。
直接法主要是基于稀疏矩阵分解技术或消去技术，采用前代回代方法进行求解。并且，由于直接法具有计算精度高、算法通用性强的优势，成为目前求解大规模稀疏线性方程组的主要手段。Davis等[30]和Yeralan等[31]概述了求解线性方程组的直接方法，包括LU、QR、Cholesky和其他因子分解、前代和回代以及相关矩阵运算。George等[32]提出了一种GPU上对称正定矩阵的多因子分解[33]方法；Rennich等[34]提出了一种超节点的Cholesky分解算法，该算法通过GPU传输消除树[35]的分支(以叶子结尾的子树)，并完全在GPU上对每个分支进行分解。Davis等[36]提出了一种基于GPU的多额稀疏QR因子分解方法[37]；Buttari等[38]对多额稀疏QR因子分解方法进行了改进，提出了一种细粒度的多额稀疏QR因子分解方法，但是细粒度的划分会影响计算的速度。
2 基于变量相关性分解的并行求解算法构造7×7稀疏线性方程组：
| $\left[\begin{array}{ccccc}A_{11} & & & A_{15} & & A_{17} \\& A_{22} & A_{23} & A_{24} & & & \\& A_{32} & A_{33} & & & A_{36} & \\& A_{12} & & A_{44} & A_{45} & & A_{47} \\A_{51} & & & A_{54} & A_{55} & & \\& & A_{63} & & & A_{66} & \\A_{71} & & & A_{74} & & A_{77}\end{array}\right] \times\left[\begin{array}{c}x_{1} \\x_{2} \\x_{3} \\x_{4} \\x_{5} \\x_{6} \\x_{7}\end{array}\right]=\left[\begin{array}{c}b_{1} \\b_{2} \\b_{3} \\b_{4} \\b_{5} \\b_{6} \\b_{7}\end{array}\right].$ |
前代过程：
| $\left[\begin{array}{ccccc}L_{11} & & & & \\& L_{22} & & & \\& L_{32} & L_{33} & & \\& L_{42} & & L_{44} & & \\L_{51} & & & L_{54} & L_{55} & \\& & L_{63} & & L_{66} & \\L_{71} & & & L_{74} & & y_{77}\end{array}\right]=\left[\begin{array}{c}y_{1} \\y_{2} \\y_{3} \\y_{4} \\y_{5} \\y_{6} \\y_{7}\end{array}\right]\left[\begin{array}{c}b_{1} \\b_{2} \\b_{3} \\b_{4} \\b_{5} \\b_{6} \\b_{7}\end{array}\right].$ |
| $\left[\begin{array}{cccc}U_{11} & & U_{15} & & U_{17} \\& U_{22} & U_{23} & U_{24} & & & \\& U_{33} & & & U_{36} & \\& & U_{44} & U_{45} & & U_{47} \\& & & U_{55} & & \\& & & U_{66} & \\& & & U_{77}\end{array}\right] \times\left[\begin{array}{c}x_{1} \\x_{2} \\x_{3} \\x_{4} \\x_{5} \\x_{6} \\x_{7}\end{array}\right]=\left[\begin{array}{c}y_{1} \\y_{2} \\y_{3} \\y_{4} \\y_{5} \\y_{6} \\y_{7}\end{array}\right].$ |
| ${y_j} = \left( {{b_j} - \sum\limits_{i = 1}^{i - 1} {{L_{ij}}} \times {y_i}} \right)/{L_{ij}}, 1 \le j \le 7.$ | (1) |
回代过程由后向前遍历，计算公式如下：
| ${x_m} = \left( {{y_m} - \sum\limits_{n = m + 1}^T {{U_{mn}}} \times {x_n}} \right)/{U_{mn}}, 1 \le m \le 7.$ | (2) |
通过上述变量求解的依赖关系，分别构造表示变量求解顺序的关联图G(V, E)，如图 1所示。其中：V代表关联图中所有节点的集合，每个节点代表一个变量；E是图中所有有向边的集合，每条有向边代表 2个相关变量求解顺序的依赖关系。有向边的起始端变量i是末端变量j的前驱变量，末端变量j是起始端变量i的后继变量。系数矩阵中第j行i列(j>i)的非零元素lij由图中节点i到节点j的有向边(i, j)表示。
 |
| 图 1 变量求解顺序的关联图G (V, E) |
| 图选项 |
以y3和x3的计算为例，验证算法的可行性和正确性。应用直接法分别求解y3和x3，
| ${{L_{32}}{y_2} + {L_{33}}{y_3} = {b_3}, }$ | (3) |
| ${{U_{33}}{x_3} + {U_{36}}{x_6} = {y_3}.}$ | (4) |
| ${{y_3} = \left( {{b_3} - {L_{32}}{y_2}} \right)/{L_{33}}, }$ | (5) |
| ${{x_3} = \left( {{y_3} - {U_{36}}{x_6}} \right)/{U_{33}}.}$ | (6) |
代入式(5)得
| $y_{3}=\left(b_{3}-L_{32} b_{2} / L_{22}\right) / L_{33}.$ | (7) |
又y6由直接法求得，y6=(b6－L63y3)/L66，将y3代入，得
| $y_{6}=\left[b_{6}-L_{63}\left(b_{3}-L_{32} b_{2} / L_{22}\right)\right] / L_{66}.$ | (8) |
| $y_{3}=b_{3}-L_{32} b_{2}, $ | (9) |
| $\begin{array}{l}x_{3}=b_{3} / U_{33}+\left[\left(-U_{36}\right) \times\left(-L_{63}\right) \times b_{3} / U_{66}\right] / U_{33}+ \\\quad\left(-L_{32} \times b_{2}\right) / U_{33}+\left[\left(-U_{36}\right) \times b_{6} / U_{66}\right] / U_{33}+ \\\quad\left[\left(-U_{36}\right) \times\left(-L_{63}\right) \times\left(-L_{32} \times b_{2}\right) / U_{66}\right] / U_{33}.\end{array}$ | (10) |
在本文算法中求解y3和x3，首先分解变量y和x的关联图，得到与节点3相连的子图如2a和2b所示。
由图 2a中y3求解子图可以看到，指向节点3的有向边只有节点2，即变量y3的前继变量只有y2，则y3的部分值由b2和b3决定。同理，图 2b中x3求解子图中，指向节点3的有向边只有一个节点6，即变量x3的前继变量只有x6，则x3的部分值由y3和y6决定。图2c的y6求解子图中，指向节点6的有向边只有节点3，即只要有y3就可计算y6，则y6的部分值由b2、b3和b6决定。综上所述，求解y3和x3只需知道b2、b3和b6。计算过程如图 3所示。
 |
| 图 2 y3、x3、y6求解子图 |
| 图选项 |
 |
| 图 3 y3、x3部分值计算过程 |
| 图选项 |
图 3a表示当且仅当b2≠0时，由Ly=b可求得y3和y6的部分值，再根据Ux=y求得x3的部分值。图 3b表示当且仅当b3≠0时，由Ly=b可求得y3和y6的部分值，再根据Ux=y求得x3的部分值。图 3c表示当且仅当b6≠0时，由Ly=b可求得y6的部分值，再根据Ux=y求得x6和x3的部分值。
由以上分析可知，y3的部分值由2部分组成，x3的部分值由5部分组成，并且所有部分值并行计算，将y3和x3的所有部分值相加得到最终值，经验证与直接法求得的y3和x3最终值相同。
3 并行算法的实现本文在GPU实验平台下，设置2N个线程，前N个线程计算y，后N个线程计算x。将求解稀疏线性方程组的过程划分成多个计算任务，且每个任务安排到一个线程计算。具体到每个线程，首先计算变量y的部分值，同时，通过y和x的相关性将上述计算的y的部分值传到变量x的相应列的线程内计算x的部分值，最后将所有部分值相加得到最终值。
由于y和x的部分值求解过程相同，本文以y的部分值求解为例说明计算过程。首先计算下三角矩阵L每一列的第1个变量的部分值，再继续计算这一列和第1个变量直接相关的变量部分值，并传到存储相应的变量部分值的线程内，再由变量部分值计算其余相关变量的部分值，以上计算过程并行执行。为了不丢失变量的部分值，引入变量部分值的计算次数。首先，线程计算与起始变量直接相关的变量，得到后继变量的部分值，同时后继变量部分值的计算次数加1。然后，线程通过循环判断某一列中起始变量部分值的计算次数，决定是否对该列进行计算。若起始变量部分值的计算次数为0，代表该列中后继变量的部分值计算完成，线程不再对该列计算；若非0，则线程继续计算该列中后继变量的部分值，并且后继变量部分值的计算次数继续加1，起始变量部分值的计算次数减1。直到计算次数减为0，不再计算并得到变量的最终值。
在并行计算过程中，多个并行线程有可能对一个变量同时进行部分值的读写或计算。所以，在应用CUDA C编程实现本文并行求解算法时，使用atomic原子操作对变量部分值的计算、累加和计算次数的递减和递增进行处理。atomic原子操作是对一个变量进行“读取-计算-写入”这3项操作中一个最小单位的执行过程。当某一个线程对变量执行原子操作时，不允许其他线程对该变量进行读写或计算的操作。基于上述机制，atomic原子操作实现了多个线程间对变量的互斥保护，确保变量操作的正确性。并且，原子操作是硬件机制的同步方式，在软件层面无需同步。设b2≠0，利用atomic原子操作计算部分值过程如图 4所示。
 |
| 图 4 b2≠0的部分值计算 |
| 图选项 |
本文以求解y为例，内核函数伪代码如图 5所示。输入数据为按列压缩存储格式的稀疏矩阵和右端项B，并设置3个n维数组y、party、ytime，分别存储y的累加值、y的部分值、y的计算次数。
 |
| 图 5 并行求解y的伪代码 |
| 图选项 |
其中，3—12行是计算和起始变量直接相关的变量的部分值过程，13—24行是由已有变量的部分值继续计算部分值的过程，上述计算过程并行执行，25行为y的最终值。
4 实验结果与分析实验平台如下：GPU为GeForce GTX 1080 Ti(Pascal)，在Ubuntu16.04 64位操作系统下，使用CUDA[39] 9.1.85版本，应用CUDA C[40]语言编程实现本文算法。
本文选取99个来自不同领域的Florida矩阵的稀疏矩阵进行测试，本文算法计算时间和cusparse库函数计算时间如图 6所示，本文算法速度提升百分比如图 6所示。图 6数据的计算时间不包含系数矩阵A和右端项B的输入时间以及计算变量的输出时间，只是对求解过程的计算时间进行比较。其中调用cusparse库函数计算时，存在分析和求解2个阶段，它的计算时间是分析时间与求解时间之和；本文算法的计算时间是变量部分值的计算时间与部分值相加时间之和。横坐标代表矩阵的行/列，纵坐标代表计算时间，单位ms。散点图直观、清晰地显示了本文算法的计算时间和cusparse库函数的计算时间。
 |
| 图 6 (网络版彩图)本文算法和cusparse库函数计算时间 |
| 图选项 |
由图 6可以清晰地看到，本文算法的计算时间明显要比cusparse库函数的计算时间短。主要原因是本文算法没有分析阶段，可以直接对变量进行求解无需等待，而cusparse库函数求解时分析时间较长，而且分析时间占总计算时间比例大，进而导致总计算时间较长，如表 1所示。
表 1 求解时间和分析时间占总时间百分比
| 矩阵规模 | 分析时间/ms | 求解时间/ms | 求解时间占比/% | 分析时间占比/% |
| 4 147 110 | 58 850 | 81.504 9 | 0.13 | 99.8 |
| 12 057 441 | 23 074.4 | 19.900 6 | 0.08 | 99.9 |
| 5 154 859 | 22 715.8 | 15.169 7 | 0.06 | 99.9 |
| 23 947 347 | 41 921.5 | 32.372 2 | 0.07 | 99.9 |
| 16 777 216 | 41 309.2 | 44.787 8 | 0.10 | 99.8 |
| 11 950 757 | 21 806.5 | 45.331 7 | 0.2 | 99.7 |
| 1 564 794 | 21 087.9 | 46.121 7 | 0.21 | 99.7 |
| 9 845 725 | 19 337.3 | 30.241 6 | 0.15 | 99.8 |
| 1 470 152 | 15 923.6 | 42.832 5 | 0.26 | 99.7 |
| 5 558 326 | 15 228.5 | 11.050 8 | 0.07 | 99.9 |
| 6 815 744 | 14 257.6 | 24.639 9 | 0.17 | 99.8 |
| 943 695 | 13 996.3 | 29.744 3 | 0.21 | 99.7 |
| 1 391 349 | 12 197.5 | 27.591 3 | 0.22 | 99.7 |
| 1 437 960 | 11 956.9 | 29.801 8 | 0.24 | 99.7 |
| 1 498 023 | 11 761.3 | 77.057 8 | 0.65 | 99.3 |
| 726 713 | 8 234.35 | 2 072.200 | 20.10 | 79.80 |
| 3 799 275 | 9 837.28 | 77.844 1 | 0.78 | 99.2 |
表选项
与cusparse库函数求解方法相比，本文算法速度提升百分比如图 7所示。加速比的计算公式为
| $\alpha=\frac{T_{\mathrm{cus}}-T}{T_{\mathrm{cus}}} \times 100 \%.$ | (11) |
 |
| 图 7 本文算法速度提升百分比 |
| 图选项 |
其中，Tcus为cusparse库函数计算时间，T为本文算法计算时间。
由图 7可以看到，本文算法速度提升百分比均在50%以上，最大可以达到99.9%，充分验证了本文算法的可行性和高效性。
5 结论本文通过分解变量之间的相关性，把求解稀疏线性方程组的过程划分为多个独立且并行的计算任务，而且变量的计算无需等待其前继变量计算全部计算完成就可进行部分值的计算，提高了算法并行度和求解速度。最后，利用Florida稀疏矩阵集构造了稀疏线性方程组，并应用本文提出的分解变量相关性的方法进行求解，同时与调用cusparse库函数的求解方法进行了比较，验证了本文并行算法的可行性和高效性。
参考文献
| [1] | SALINAS S, LUO C Q, CHEN X H, et al. Efficient secure outsourcing of large-scale sparse linear systems of equations[J]. IEEE Transactions on Big Data, 2017, 4(1): 26-39. |
| [2] | 沈海龙.线性代数系统迭代解法与预条件方法研究[D].沈阳: 东北大学, 2013. SHEN H L. The study of iterative methods and preconditioning methods for linear algebraic systems[D]. Shenyang: Northeastern University, 2013. (in Chinese) |
| [3] | 周挺辉, 赵文恺, 严正, 等. 基于图形处理器的电力系统稀疏线性方程组求解方法[J]. 电力系统自动化, 2015, 39(2): 74-80. ZHOU T K, ZHAO W K, YAN Z, et al. A method for solving sparse linear equations of power systems based on GPU[J]. Automation of Electric Power System, 2015, 39(2): 74-80. (in Chinese) |
| [4] | 谢力, 王武, 冯仰德. 基于多层半可分结构矩阵的快速算法与并行实现[J]. 数值计算与计算机应用, 2017, 38(1): 37-48. XIE L, WANG W, FENG Y D. A fast algorithm based on hierarchically semiseparable structured matrix and its parallel implementation[J]. Journal on Numerical Methods and Computer Applications, 2017, 38(1): 37-48. (in Chinese) |
| [5] | KECKLER S W, DALLY W J, KHAILANY B, et al. GPUs and the future of parallel computing[J]. IEEE Micro, 2011, 31(5): 7-17. DOI:10.1109/MM.2011.89 |
| [6] | 郑方, 张昆, 邬贵明, 等. 面向高性能计算的众核处理器结构级高能效技术[J]. 计算机学报, 2014, 37(10): 2176-2186. ZHENG F, ZHANG K, WU G M, et al. Architecture techniques of many-core processor for energy-efficient in high performance computing[J]. Journal of Computers, 2014, 37(10): 2176-2186. (in Chinese) |
| [7] | 朱宇兰. 基于GPU通用计算的并行算法和计算框架的实现[J]. 山东农业大学学报(自然科学版), 2016, 47(3): 473-476, 480. ZHU Y L. Parallel algorithm based on general purpose computing on GPU and the implementation of calculation framework[J]. Journal of Shandong Agricultural University (Natural Science Edition), 2016, 47(3): 473-476, 480. DOI:10.3969/j.issn.1000-2324.2016.03.030 (in Chinese) |
| [8] | PARK J, SMELYANSKIY M, SUNDARAM N, et al. Sparsifying synchronization for high-performance shared memory sparse-triangular solver[M]//KUNKEL J M, LUDWIG T, MEUER H W. Supercomputing. Cham: Springer, 2014: 124-140. |
| [9] | NAUMOV M. Parallel solution of sparse triangular linear systems in the preconditioned iterative methods on the GPU[R]. Technical Report NVR-2011-001 NVIDIA. NVIDIA, 2011. |
| [10] | LI R P, SAAD Y. GPU-accelerated preconditioned iterative linear solvers[J]. The Journal of Supercomputing, 2013, 63(2): 443-466. DOI:10.1007/s11227-012-0825-3 |
| [11] | WOLF M M, HEROUX M A, BOMAN E G. Factors impacting performance of multithreaded sparse triangular solve[M]//PALMA J M L M, DAYDé M, MARQUES O, et al. High Performance Computing for Computational Science-VECPAR 2010. Berlin, Heidelberg: Springer, 2011, 6449: 32-44. |
| [12] | KABIR H, BOOTH J D, AUPY G, et al. STS-k: A multilevel sparse triangular solution scheme for NUMA multicores[C]//Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. New York, USA: ACM, 2015. |
| [13] | MAADASAMY S, DORAISWAMY H, NATARAJAN V. A hybrid parallel algorithm for computing and tracking level set topology[C]//2012 19th International Conference on High Performance Computing. Pune, India: IEEE, 2012. |
| [14] | PICCIAU A, INGGS G E, WICKERSON J, et al. Balancing locality and concurrency: Solving sparse triangular systems on GPUs[C]//2016 IEEE 23rd International Conference on High Performance Computing. Hyderabad, India: IEEE, 2016. |
| [15] | SUCHOSKI B, SEVERN C, SHANTHARAM M, et al. Adapting sparse triangular solution to GPUs[C]//201241stInternational Conferenceon Parallel ProcessingWorkshops. Pittsburgh, USA: IEEE, 2012. |
| [16] | LI A, VAN DEN BRAAK G J, CORPORAAL H, et al. Fine-grained synchronizations and dataflow programming on GPUs[C]//Proceedings of the 29th ACM on International Conference on Supercomputing. New York, USA: ACM, 2015. |
| [17] | CHE S, SHEAFFER J W, SKADRON K. Dymaxion: optimizing memory access patterns for heterogeneous systems[C]//Proceedings of the 2011 International Conference for High Performance Computing, Networking, Storage and Analysis. New York, USA: IEEE, 2011. |
| [18] | DAVIS T A, HU Y F. The university of Florida sparse matrix collection[J]. ACM Transactions on Mathematical Software, 2011, 38(1): 1-25. |
| [19] | NVIDIA Corporation. CUDA Toolkit 10. 0, cusparse library, 2018[EB/OL]. (2018-09-19) http://developer.nvidia.com/cuda-toolkit-40. |
| [20] | Ahamed A K C, MAGOULōS F. Parallel sub-structuring methods for solving sparse linear systems on a cluster of GPUs[C]//2014 IEEE International Conference on High Performance Computing and Communications, 2014 IEEE 6th International Symposium on Cyberspace Safety and Security, 2014 IEEE 11th International Conference on Embedded Software and Syst. Paris, France: IEEE, 2015. |
| [21] | SAAD Y. Iterative methods for sparse linear systems[M]. Boston: PWS Pub. Co., 2009: 216-220. |
| [22] | YOUNG D M. Iterative solution of large linear systems[M]. Amsterdam: Elsevier, 2014. |
| [23] | DAVIS T A. Direct methods for sparse linear systems[M]. Philadelphia: Fundamentals of Algorithms, 2006. |
| [24] | DAVIS T A, RAJAMANICKAM S, SID-LAKHDAR W M. A survey of direct methods for sparse linear systems[J]. Acta Numerica, 2016, 25: 383-566. DOI:10.1017/S0962492916000076 |
| [25] | AMENT M, KNITTEL G, WEISKOPF D, et al. A parallel preconditioned conjugate gradient solver for the Poisson problem on a multi-GPU platform[C]//2010 18th Euromicro Conference on Parallel, Distributed and Network-based Processing. Pisa, Italy: IEEE, 2010: 583-592. |
| [26] | BAHI J M, RAPHA?L C, KHODJA L Z. Parallel sparse linear solver GMRES for GPU clusters with compression of exchanged data[M]. Berlin Heidelberg: Parallel Processing Workshops, 2012. |
| [27] | KLIE H M, SUDAN H H, LI R, et al. Exploiting capabilities of many core platforms in reservoir simulation[C]//Society of Petroleum Engineers SPE Reservoir Simulation Symposium. The Woodlands, Texas, USA: SPE Reservoir Simulation Symposium, 2011. |
| [28] | DZIEKONSKI A, LAMECKI A, MROZOWSKI M. Jacobi and Gauss-Seidel preconditioned complex conjugate gradient method with GPU acceleration for finite element method[C]//The 40th European Microwave Conference. Paris, France: IEEE, 2010. |
| [29] | GAIKWAD A, TOKE I M. Parallel iterative linear solvers on GPU: A financial engineering case[C]//2010 18th Euromicro Conference on Parallel, Distributed and Network-based Processing. Pisa, Italy: IEEE, 2010. |
| [30] | DUFF I S, ERISMAN A M, REID J K. Direct methods for sparse matrices[M]. 2nd ed. Oxford: Oxford University Press, 2017. |
| [31] | Yeralan S N, Davis T A, Ranka S. Sparse QR factorization on GPU architectures[R]. Florida: University of Florida, 2013: 1-36. |
| [32] | AGULLO E, BUTTARI A, GUERMOUCHE A, et al. Multifrontal QR factorization for multicore architectures over runtime systems[M]//WOLF F, MOHR B, AN MEY D. Euro-Par 2013 Parallel Processing. Berlin Heidelberg: Springer, 2013. |
| [33] | BUTTARI A. Fine-grained multithreading for themultifrontal QR factorization of sparse matrices[J]. SIAM Journal on Scientific Computing, 2013, 35(4). |
| [34] | GEORGE T, SAXENAA V, GUPTA A, et al. Multifrontal factorization of sparse SPD matrices on GPUs[C]//2011 IEEE International Parallel & Distributed Processing Symposium, Anchorage, USA: IEEE, 2011: 372-383. |
| [35] | LUCAS R F, WAGENBRETH G, TRAN J J, et al. Multifrontal sparse matrix factorization on graphics processing units[EB/OL].[2012-01-13]. http://dx.doi.org/.2012. |
| [36] | RENNICH S C, STOSIC D, DAVIS T A. Accelerating sparse cholesky factorization on GPUs[J]. Parallel Computing, 2016, 59: 140-150. DOI:10.1016/j.parco.2016.06.004 |
| [37] | CHEN X M, REN L, WANG Y, et al. GPU-accelerated sparse LU factorization for circuit simulation with performance modeling[J]. IEEE Transactions on Parallel and Distributed Systems, 2015, 26(3): 786-795. DOI:10.1109/TPDS.2014.2312199 |
| [38] | KAYAASLAN E, U?AR B. Reducing elimination tree height for parallel LU factorization of sparse unsymmetric matrices[C]//201421st International Conference on High Performance Computing. Dona Paula, India: IEEE, 2014: 1-10. |
| [39] | NVIDIA. CUDA toolkit[EB/OL].[2019-05-18]. https://developer.nvidia.com/cuda-toolkit. |
| [40] | NVIDIA. NVIDIA CUDA C programming guide. Version 3.1[M]. San Jose: NVIDIA, 2010. |
