使用cuSPARSE和cuBLAS的不完全LU分解与Cholesky分解预处理迭代方法

白皮书描述了如何利用cuSPARSE和cuBLAS库在不完全-LU和Cholesky预处理迭代方法中实现比CPU快2倍的速度。

1. 简介

大规模稀疏线性系统的求解是计算力学、大气建模、地球物理学、生物学、电路仿真以及计算科学与工程领域许多其他应用中的重要问题。通常，这些线性系统可以使用直接法或预处理迭代法来求解。尽管直接法通常更可靠，但它们通常需要较大的内存需求，并且在大规模并行计算机平台上的扩展性不佳。

迭代方法更适合并行化，因此可用于解决更大规模的问题。目前最流行的迭代方案属于Krylov子空间方法家族，包括分别用于非对称和对称正定(s.p.d.)线性系统的双共轭梯度稳定法(BiCGStab)和共轭梯度法(CG)[2], [11]。我们将在下一节更详细地描述这些方法。

在实践中，我们经常使用各种预处理技术来提高迭代方法的收敛性。本白皮书重点讨论不完全-LU和Cholesky预处理[11]，这是最流行的预处理技术之一。它计算系数矩阵的不完全分解，并在迭代方法的每次迭代中需要求解下三角和上三角线性方程组。

为了实现预条件化的BiCGStab和CG算法，我们使用了cuSPARSE库中实现的稀疏矩阵-向量乘法[3]、[15]以及稀疏三角求解[8]、[16]。需要指出的是，这些算法的底层实现充分利用了CUDA并行编程范式[5]、[9]、[13]的优势，使我们能够充分挖掘图形处理器(GPU)的计算资源。在我们的数值实验中，不完全分解是在CPU（主机端）上执行的，生成的下三角和上三角因子会在开始迭代方法前传输到GPU（设备端）内存。不过，不完全分解的计算过程同样可以在GPU上获得加速。

我们指出，这些迭代方法中可用的并行性高度取决于当前系数矩阵的稀疏模式。在我们的数值实验中，使用GPU上的cuSPARSE和cuBLAS库进行不完全-LU和Cholesky预处理的迭代方法，相比CPU上的MKL[17]实现平均获得了超过2倍的加速。例如，图1展示了来自不同应用场景的矩阵在使用0填充(ilu0)的不完全-LU和Cholesky分解预处理迭代方法时的加速效果。最后一节将对此进行更详细的说明。

不完全LU Cholesky（零填充）预处理迭代方法的加速比

在接下来的章节中，我们将简要介绍相关方法，并讨论并行稀疏矩阵-向量乘法和三角求解算法在这些方法中所起的作用。

2. 预处理迭代方法

让我们考虑这个线性系统

\(A\mathbf{x} = \mathbf{f}\)

其中 \(A \in \mathbb{R}^{n \times n}\) 是一个非奇异系数矩阵，而 \(\mathbf{x},\mathbf{f} \in \mathbb{R}^{n}\) 分别是解向量和右端向量。

一般来说，迭代方法从一个初始猜测开始，通过一系列步骤逐步逼近更精确的解。迭代方法主要分为两类：(i) 静态迭代方法，如基于分裂的Jacobi和Gauss-Seidel(GS)；(ii) 非静态迭代方法，如Krylov子空间方法家族，包括CG和BiCGStab。如前所述，本白皮书将重点讨论后者。

迭代方法的收敛性高度依赖于系数矩阵的谱特性，通过使用预处理可以显著改善收敛效果。预处理通过修改线性系统系数矩阵的谱特性，以减少收敛所需的迭代步数。通常需要找到一个预处理矩阵\(M\)，使得\(M^{- 1}\)能较好地近似\(A^{- 1}\)，并且用\(M\)求解系统相对容易。

对于对称正定矩阵\(A\)，我们可以令\(M\)为其不完全Cholesky分解，使得\(A \approx M = {\widetilde{R}}^{T}\widetilde{R}\)，其中\(\widetilde{R}\)是上三角矩阵。假设\(M\)非奇异，则\({\widetilde{R}}^{- T}A{\widetilde{R}}^{- 1}\)也是对称正定的，此时我们可以求解预处理后的线性方程组来代替原线性方程组(1)的求解。

\(\left( {{\widetilde{R}}^{- T}A{\widetilde{R}}^{- 1}} \right)\left( {\widetilde{R}\mathbf{x}} \right) = {\widetilde{R}}^{- T}\mathbf{f}\)

预处理共轭梯度迭代法的伪代码如算法1所示。

2.1. 算法1 共轭梯度法(CG)

1:	\(\text{设初始猜测为 }\mathbf{x}_{0}\text{, 计算 }\mathbf{r}\leftarrow\mathbf{f} - A\mathbf{x}_{0}\)
2:	\(\textbf{for }i\leftarrow 1,2,...\text{ until convergence }\textbf{do}\)
3:	\(\quad\quad\text{Solve }M\mathbf{z}\leftarrow\mathbf{r}\)	\(\vartriangleright \text{稀疏下三角和上三角求解}\)
4:	\(\quad\quad\rho_{i}\leftarrow\mathbf{r}^{T}\mathbf{z}\)
5:	\(\quad\quad\textbf{if }i==1\textbf{ then}\)
6:	\(\quad\quad\quad\quad\mathbf{p}\leftarrow\mathbf{z}\)
7:	\(\quad\quad\textbf{否则}\)
8:	\(\quad\quad\quad\quad\beta\leftarrow\frac{\rho_{i}}{\rho_{i - 1}}\)
9:	\(\quad\quad\quad\quad\mathbf{p}\leftarrow\mathbf{z} + \beta\mathbf{p}\)
10:	\(\quad\quad\textbf{end if}\)
11:	\(\quad\quad\text{计算 }\mathbf{q}\leftarrow A\mathbf{p}\)	\(\vartriangleright \text{稀疏矩阵-向量乘法}\)
12:	\(\quad\quad\alpha\leftarrow\frac{\rho_{i}}{\mathbf{p}^{T}\mathbf{q}}\)
13:	\(\quad\quad\mathbf{x}\leftarrow\mathbf{x} + \alpha\mathbf{p}\)
14:	\(\quad\quad\mathbf{r}\leftarrow\mathbf{r} - \alpha\mathbf{q}\)
15:	\(\textbf{end for}\)

请注意，在不完全Cholesky预处理共轭梯度(CG)迭代方法的每次迭代中，我们需要执行一次稀疏矩阵向量乘法和两次三角求解。下面展示了使用C编程语言中cuSPARSE和cuBLAS库的相应CG代码。

/***** CG Code *****/
/* ASSUMPTIONS:
   1. The cuSPARSE and cuBLAS libraries have been initialized.
   2. The appropriate memory has been allocated and set to zero.
   3. The matrix A (valA, csrRowPtrA, csrColIndA) and the incomplete-
      Cholesky upper triangular factor R (valR, csrRowPtrR, csrColIndR)
      have been computed and are present in the device (GPU) memory. */

//create the info and analyse the lower and upper triangular factors
cusparseCreateSolveAnalysisInfo(&inforRt);
cusparseCreateSolveAnalysisInfo(&inforR);
cusparseDcsrsv_analysis(handle,CUSPARSE_OPERATION_TRANSPOSE,
                      n, descrR, valR, csrRowPtrR, csrColIndR, inforRt);
cusparseDcsrsv_analysis(handle,CUSPARSE_OPERATION_NON_TRANSPOSE,
                      n, descrR, valR, csrRowPtrR, csrColIndR, inforR);

//1: compute initial residual r = f -  A x0 (using initial guess in x)
cusparseDcsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, n, n, 1.0,
               descrA, valA, csrRowPtrA, csrColIndA, x, 0.0, r);
cublasDscal(n,-1.0, r, 1);
cublasDaxpy(n, 1.0, f, 1, r, 1);
nrmr0 = cublasDnrm2(n, r, 1);

//2: repeat until convergence (based on max. it. and relative residual)
for (i=0; i<maxit; i++){
    //3: Solve M z = r (sparse lower and upper triangular solves)
    cusparseDcsrsv_solve(handle, CUSPARSE_OPERATION_TRANSPOSE,
                         n, 1.0, descrpR, valR, csrRowPtrR, csrColIndR,
                         inforRt, r, t);
    cusparseDcsrsv_solve(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
                         n, 1.0, descrpR, valR, csrRowPtrR, csrColIndR,
                         inforR, t, z);

    //4: \rho = r^{T} z
    rhop= rho;
    rho = cublasDdot(n, r, 1, z, 1);
    if (i == 0){
        //6: p = z
        cublasDcopy(n, z, 1, p, 1);
    }
    else{
        //8: \beta = rho_{i} / \rho_{i-1}
        beta= rho/rhop;
        //9: p = z + \beta p
        cublasDaxpy(n, beta, p, 1, z, 1);
        cublasDcopy(n, z, 1, p, 1);
    }

    //11: Compute q = A p (sparse matrix-vector multiplication)
    cusparseDcsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, n, n, 1.0,
                   descrA, valA, csrRowPtrA, csrColIndA, p, 0.0, q);

    //12: \alpha = \rho_{i} / (p^{T} q)
    temp = cublasDdot(n, p, 1, q, 1);
    alpha= rho/temp;
    //13: x = x + \alpha p
    cublasDaxpy(n, alpha, p, 1, x, 1);
    //14: r = r - \alpha q
    cublasDaxpy(n,-alpha, q, 1, r, 1);

    //check for convergence
    nrmr = cublasDnrm2(n, r, 1);
    if (nrmr/nrmr0 < tol){
        break;
    }
}

//destroy the analysis info (for lower and upper triangular factors)
cusparseDestroySolveAnalysisInfo(inforRt);
cusparseDestroySolveAnalysisInfo(inforR);

对于非对称矩阵\(A\)，我们可以令\(M\)为其不完全LU分解，使得\(A \nvdash M = \widetilde{L}\widetilde{U}\)，其中\(\widetilde{L}\)和\(\widetilde{U}\)分别为下三角和上三角矩阵。假设\(M\)是非奇异的，那么\(M^{- 1}A\)也是非奇异的，我们可以求解预处理后的线性系统来代替求解原线性系统(1)。

\(\left( {M^{- 1}A} \right)\mathbf{x} = M^{- 1}\mathbf{f}\)

预处理BiCGStab迭代方法的伪代码如算法2所示。

2.2. 算法2 双共轭梯度稳定法(BiCGStab)

1:	\(\text{设初始猜测为 }\mathbf{x}_{0}\text{, 计算 }\mathbf{r}\leftarrow\mathbf{f} - A\mathbf{x}_{0}\)
2:	\(\text{设置 }\mathbf{p}\leftarrow\mathbf{r}\text{ 并选择 }\widetilde{\mathbf{r}}\text{，例如可以设置 }\widetilde{\mathbf{r}}\leftarrow\mathbf{r}\)
3:	\(\textbf{for }i\leftarrow 1,2,...\text{ until convergence }\textbf{do}\)
4:	\(\quad\quad\rho_{i}\leftarrow{\widetilde{\mathbf{r}}}^{T}\mathbf{r}\)
5:	\(\quad\quad\textbf{if }\rho_{i}==0.0\textbf{ then}\)
6:	\(\quad\quad\quad\quad\text{方法失败}\)
7:	\(\quad\quad\textbf{end if}\)
8:	\(\quad\quad\textbf{if }i > 1\textbf{ then}\)
9:	\(\quad\quad\quad\quad\textbf{if }\omega==0.0\textbf{ then}\)
10:	\(\quad\quad\quad\quad\quad\quad\text{方法失败}\)
11:	\(\quad\quad\textbf{end if}\)
12:	\(\quad\quad\quad\quad\beta\leftarrow\frac{\rho_{i}}{\rho_{i - 1}} times \frac{\alpha}{\omega}\)
13:	\(\quad\quad\quad\quad\mathbf{p}\leftarrow\mathbf{r} + \beta\left( {\mathbf{p} - \omega\mathbf{v}} \right)\)
14:	\(\quad\quad\textbf{end if}\)
15:	\(\quad\quad\text{Solve }M\hat{\mathbf{p}}\leftarrow\mathbf{p}\)	\(\vartriangleright \text{稀疏上下三角矩阵求解}\)
16:	\(\quad\quad\text{计算 }\mathbf{q}\leftarrow A\hat{\mathbf{p}}\)	\(\vartriangleright \text{稀疏矩阵-向量乘法}\)
17:	\(\quad\quad\alpha\leftarrow\frac{\rho_{i}}{{\widetilde{\mathbf{r}}}^{T}\mathbf{q}}\)
18:	\(\quad\quad\mathbf{s}\leftarrow\mathbf{r} - \alpha\mathbf{q}\)
19:	\(\quad\quad\mathbf{x}\leftarrow\mathbf{x} + \alpha\hat{\mathbf{p}}\)
20:	\(\quad\quad\textbf{if }\left\| s \right\|_{2} \leq \mathit{tol}\textbf{ then}\)
21:	\(\quad\quad\quad\quad\text{方法已收敛}\)
22:	\(\quad\quad\textbf{end if}\)
23:	\(\quad\quad\text{求解 }M\hat{\mathbf{s}}\leftarrow\mathbf{s}\)	\(\vartriangleright \text{稀疏上下三角矩阵求解}\)
24:	\(\quad\quad\text{计算 }\mathbf{t}\leftarrow A\hat{\mathbf{s}}\)	\(\vartriangleright \text{稀疏矩阵-向量乘法}\)
25:	\(\quad\quad\omega\leftarrow\frac{\mathbf{t}^{T}\mathbf{s}}{\mathbf{t}^{T}\mathbf{t}}\)
26:	\(\quad\quad\mathbf{x}\leftarrow\mathbf{x} + \omega\hat{\mathbf{s}}\)
27:	\(\quad\quad\mathbf{r}\leftarrow\mathbf{s} - \omega\mathbf{t}\)
28:	\(\textbf{end for}\)

请注意，在不完全LU预处理的BiCGStab迭代方法的每次迭代中，我们需要执行两次稀疏矩阵向量乘法和四次三角求解。下面展示了使用C编程语言中cuSPARSE和cuBLAS库的相应BiCGStab代码。

/***** BiCGStab Code *****/
/* ASSUMPTIONS:
   1. The cuSPARSE and cuBLAS libraries have been initialized.
   2. The appropriate memory has been allocated and set to zero.
   3. The matrix A (valA, csrRowPtrA, csrColIndA) and the incomplete-
      LU lower L (valL, csrRowPtrL, csrColIndL)  and upper U (valU,
      csrRowPtrU, csrColIndU) triangular factors have been
      computed and are present in the device (GPU) memory. */

//create the info and analyse the lower and upper triangular factors
cusparseCreateSolveAnalysisInfo(&infoL);
cusparseCreateSolveAnalysisInfo(&infoU);
cusparseDcsrsv_analysis(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
                        n, descrL, valL, csrRowPtrL, csrColIndL, infoL);
cusparseDcsrsv_analysis(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
                        n, descrU, valU, csrRowPtrU, csrColIndU, infoU);

//1: compute initial residual r = b - A x0 (using initial guess in x)
cusparseDcsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, n, n, 1.0,
               descrA, valA, csrRowPtrA, csrColIndA, x, 0.0, r);
cublasDscal(n,-1.0, r, 1);
cublasDaxpy(n, 1.0, f, 1, r, 1);
//2: Set p=r and \tilde{r}=r
cublasDcopy(n, r, 1, p, 1);
cublasDcopy(n, r, 1, rw,1);
nrmr0 = cublasDnrm2(n, r, 1);

//3: repeat until convergence (based on max. it. and relative residual)
for (i=0; i<maxit; i++){
    //4: \rho = \tilde{r}^{T} r
    rhop= rho;
    rho = cublasDdot(n, rw, 1, r, 1);
    if (i > 0){
        //12: \beta = (\rho_{i} / \rho_{i-1}) ( \alpha / \omega )
        beta= (rho/rhop)*(alpha/omega);
        //13: p = r + \beta (p - \omega v)
        cublasDaxpy(n,-omega,q, 1, p, 1);
        cublasDscal(n, beta, p, 1);
        cublasDaxpy(n, 1.0,  r, 1, p, 1);
    }
    //15: M \hat{p} = p (sparse lower and upper triangular solves)
    cusparseDcsrsv_solve(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
                         n, 1.0, descrL, valL, csrRowPtrL, csrColIndL,
                         infoL, p, t);
    cusparseDcsrsv_solve(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
                         n, 1.0, descrU, valU, csrRowPtrU, csrColIndU,
                         infoU, t, ph);

    //16: q = A \hat{p} (sparse matrix-vector multiplication)
    cusparseDcsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, n, n, 1.0,
                   descrA, valA, csrRowPtrA, csrColIndA, ph, 0.0, q);

    //17: \alpha = \rho_{i} / (\tilde{r}^{T} q)
    temp = cublasDdot(n, rw, 1, q, 1);
    alpha= rho/temp;
    //18: s = r - \alpha q
    cublasDaxpy(n,-alpha, q, 1, r, 1);
    //19: x = x + \alpha \hat{p}
    cublasDaxpy(n, alpha, ph,1, x, 1);

    //20: check for convergence
    nrmr = cublasDnrm2(n, r, 1);
    if (nrmr/nrmr0 < tol){
        break;
    }

    //23: M \hat{s} = r (sparse lower and upper triangular solves)
    cusparseDcsrsv_solve(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
                         n, 1.0, descrL, valL, csrRowPtrL, csrColIndL,
                         infoL, r, t);
    cusparseDcsrsv_solve(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
                         n, 1.0, descrU, valU, csrRowPtrU, csrColIndU,
                         infoU, t, s);

    //24: t = A \hat{s} (sparse matrix-vector multiplication)
    cusparseDcsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, n, n, 1.0,
                   descrA, valA, csrRowPtrA, csrColIndA, s, 0.0, t);

    //25: \omega = (t^{T} s) / (t^{T} t)
    temp = cublasDdot(n, t, 1, r, 1);
    temp2= cublasDdot(n, t, 1, t, 1);
    omega= temp/temp2;
    //26: x = x + \omega \hat{s}
    cublasDaxpy(n, omega, s, 1, x, 1);
    //27: r = s - \omega t
    cublasDaxpy(n,-omega, t, 1, r, 1);

    //check for convergence
    nrmr = cublasDnrm2(n, r, 1);
    if (nrmr/nrmr0 < tol){
        break;
    }
}

//destroy the analysis info (for lower and upper triangular factors)
cusparseDestroySolveAnalysisInfo(infoL);
cusparseDestroySolveAnalysisInfo(infoU);

如图2所示，在不完全-LU和Cholesky预处理迭代方法的每次迭代中，大部分时间都花费在稀疏矩阵-向量乘法和三角求解上。稀疏矩阵-向量乘法已在以下参考文献[3]、[15]中得到广泛研究。稀疏三角求解则不太为人所知，因此我们简要指出其中用于探索并行性的策略，并建议读者参阅NVIDIA技术报告[8]获取更多细节。

预处理迭代方法在GPU上消耗的总时间拆分

要理解稀疏三角求解背后的核心思想，需注意到虽然前向和后向替换对于稠密三角系统本质上是串行算法，但对于稀疏三角系统，并不一定存在对先前获得的解元素的依赖关系。我们采用的策略正是利用这些依赖关系的缺失，将求解过程分为两个阶段，如[1]、[4]、[6]、[7]、[8]、[10]、[12]、[14]中所述。

分析阶段构建数据依赖图，根据矩阵稀疏模式将独立行分组为不同层级。求解阶段依次迭代已构建的各个层级，并行计算单个层级内所有行对应的解元素。请注意，根据设计每个层级内的行彼此独立，但都依赖于前一层级的至少一行。

分析阶段只需执行一次，通常比可以多次执行的求解阶段慢得多。这种安排特别适合不完全-LU和Cholesky预处理的迭代方法。

3. 数值实验

本节我们研究了不完全LU分解和Cholesky预处理的BiCGStab与CG迭代方法的性能表现。在数值实验中，我们选用了佛罗里达大学稀疏矩阵集[18]中的12个矩阵。7个对称正定矩阵和5个非对称矩阵按其行数(m)、列数(n=m)和非零元素数(nnz)分组，并按递增顺序展示在表1中。

表1. 对称正定(s.p.d.)与非对称测试矩阵

#	矩阵	行数,列数	非零元素数	对称正定	应用领域
1	离岸	259,789	4,242,673	是	地球物理学
2	af_shell3	504,855	17,562,051	是	力学
3	parabolic_fem	525,825	3,674,625	是	通用
4	apache2	715,176	4,817,870	是	机械
5	ecology2	999,999	4,995,991	是	生物学
6	thermal2	1,228,045	8,580,313	yes	热力模拟
7	G3_circuit	1,585,478	7,660,826	是	电路仿真
8	FEM_3D_thermal2	147,900	3,489,300	否	力学
9	thermomech_dK	204,316	2,846,228	否	力学
10	ASIC_320ks	321,671	1,316,08511	否	电路仿真
11	cage13	445,315	7,479,343	否	生物学
12	atmosmodd	1,270,432	8,814,880	no	大气模型

在以下实验中，我们使用的硬件系统配置为NVIDIA C2050（ECC启用）GPU和Intel Core i7 CPU 950 @ 3.07GHz，运行64位Linux操作系统Ubuntu 10.04 LTS，搭配cuSPARSE库4.0和MKL 10.2.3.029。未设置MKL_NUM_THREADS和MKL_DYNAMIC环境变量，以便MKL能自动选择最优线程数。

我们使用MKL例程csrilu0和csrilut分别计算不完全-LU分解和Cholesky分解，其中填充量为0和阈值填充。在csrilut例程中，我们允许三种不同级别的填充量，分别表示为(5,10^-3)、(10,10^-5)和(20,10^-7)。一般来说，\(\left( k,\mathit{tol} \right)\)填充量基于\(nnz/n + k\)每行允许的最大元素数量，并丢弃幅度小于\(\left| l_{ij} \middle| , \middle| u_{ij} \middle| < \mathit{tol} \times \left\| \mathbf{a}_{i}^{T} \right\|_{2} \right.\)的元素，其中\(l_{ij}\)、\(u_{ij}\)和\(\mathbf{a}_{i}^{T}\)分别是下三角矩阵\(L\)、上三角矩阵\(U\)的元素以及系数矩阵\(A\)的第i行。

我们比较了在GPU上使用cuSPARSE和cuBLAS库以及在CPU上使用MKL库实现的BiCGStab和CG迭代方法。在我们的实验中，初始猜测设为零，右侧向量\(\mathbf{f} = A\mathbf{e}\)其中\(\mathbf{e}^{T}{= (1,\ldots,1)}^{T}\)，停止条件设为最大迭代次数2000次或相对残差\(\left\| \mathbf{r}_{i} \right\|_{2}/\left\| \mathbf{r}_{0} \right\|_{2} < 10^{- 7}\)，其中\(\mathbf{r}_{i} = \mathbf{f} - A\mathbf{x}_{i}\)表示第i次迭代时的残差。

表2. csrilu0 预处理的共轭梯度法与双共轭梯度稳定法

	ilu0	CPU		GPU	加速比
#	实际时间(秒)	复制时间(秒)	求解时间(秒)	\(\frac{\left\| \mathbf{r}_{i} \right\|_{2}}{\left\| \mathbf{r}_{0} \right\|_{2}}\)	迭代次数	求解时间(秒)	\(\frac{\left\| \mathbf{r}_{i} \right\|_{2}}{\left\| \mathbf{r}_{0} \right\|_{2}}\)	迭代次数	对比ilu0
1	0.38	0.02	0.72	8.83E-08	25	1.52	8.83E-08	25	0.57
2	1.62	0.04	38.5	1.00E-07	569	33.9	9.69E-08	571	1.13
3	0.13	0.01	39.2	9.84E-08	1044	6.91	9.84E-08	1044	5.59
4	0.12	0.01	35.0	9.97E-08	713	12.8	9.97E-08	713	2.72
5	0.09	0.01	107	9.98E-08	1746	55.3	9.98E-08	1746	1.92
6	0.40	0.02	155	9.96E-08	1656	54.4	9.79E-08	1656	2.83
7	0.16	0.02	20.2	8.70E-08	183	8.61	8.22E-08	183	2.32
8	0.32	0.02	0.13	5.25E-08	4	0.52	5.25E-08	4	0.53
9	0.20	0.01	72.7	1.96E-04	2000	40.0	2.08E-04	2000	1.80
10	0.11	0.01	0.27	6.33E-08	6	0.12	6.33E-08	6	1.59
11	0.70	0.03	0.28	2.52E-08	2.5	0.15	2.52E-08	2.5	1.10
12	0.25	0.04	12.5	7.33E-08	76.5	4.30	9.69E-08	74.5	2.79

表3. csrilut(5,10^-3)预处理共轭梯度法与双共轭梯度稳定法

		ilut(5,10-3)		CPU		GPU		加速比
#	实际时间(秒)	复制时间(秒)	求解时间(秒)	\(\frac{\left\| \mathbf{r}_{i} \right\|_{2}}{\left\| \mathbf{r}_{0} \right\|_{2}}\)	# 迭代次数	求解时间(秒)	\(\frac{\left\| \mathbf{r}_{i} \right\|_{2}}{\left\| \mathbf{r}_{0} \right\|_{2}}\)	# 迭代次数	对比 ilut (5,10-3)	对比 ilu0
1	0.14	0.01	1.17	9.70E-08	32	1.82	9.70E-08	32	0.67	0.69
2	0.51	0.03	49.1	9.89E-08	748	33.6	9.89E-08	748	1.45	1.39
3	1.47	0.02	11.7	9.72E-08	216	6.93	9.72E-08	216	1.56	1.86
4	0.17	0.01	67.9	9.96E-08	1495	26.5	9.96E-08	1495	2.56	5.27
5	0.55	0.04	59.5	9.22E-08	653	71.6	9.22E-08	653	0.83	1.08
6	3.59	0.05	47.0	9.50E-08	401	90.1	9.64E-08	401	0.54	0.92
7	1.24	0.05	23.1	8.08E-08	153	24.8	8.08E-08	153	0.93	2.77
8	0.82	0.03	0.12	3.97E-08	2	1.12	3.97E-08	2	0.48	1.10
9	0.10	0.01	54.3	5.68E-04	2000	24.5	1.58E-04	2000	2.21	1.34
10	0.12	0.01	0.16	4.89E-08	4	0.08	6.45E-08	4	1.37	1.15
11	4.99	0.07	0.36	1.40E-08	2.5	0.37	1.40E-08	2.5	0.99	6.05
12	0.32	0.03	39.2	7.05E-08	278.5	10.6	8.82E-08	270.5	3.60	8.60

数值实验结果展示在表2至表5中，其中我们列出了GPU上迭代方法相比CPU的加速比（speedup）、收敛所需的迭代次数（# it.）、达到的相对残差（\(\frac{\left\| \mathbf{r}_{i} \right\|_{2}}{\left\| \mathbf{r}_{0} \right\|_{2}}\)）以及分解耗时（fact.）、线性系统迭代求解耗时（solve）、将上下三角因子传输至GPU的cudaMemcpy耗时（copy）。我们计算加速比时包含了不完全-LU分解和Cholesky分解的耗时，以及将三角因子从CPU内存传输至GPU内存的耗时。

表4. 使用csrilut(10,10^-5)预处理的共轭梯度法(CG)和稳定双共轭梯度法(BiCGStab)

		ilut(10,10-5)		CPU		GPU		加速比
#	实际时间(秒)	复制时间(秒)	求解时间(秒)	\(\frac{\left\| \mathbf{r}_{i} \right\|_{2}}{\left\| \mathbf{r}_{0} \right\|_{2}}\)	# 迭代次数	求解时间(秒)	\(\frac{\left\| \mathbf{r}_{i} \right\|_{2}}{\left\| \mathbf{r}_{0} \right\|_{2}}\)	# 迭代次数	对比 ilut (10,10-5)	对比 ilu0
1	0.15	0.01	1.06	8.79E-08	34	1.96	8.79E-08	34	0.57	0.63
2	0.52	0.03	60.0	9.86E-08	748	38.7	9.86E-08	748	1.54	1.70
3	3.89	0.03	9.02	9.79E-08	147	5.42	9.78E-08	147	1.38	1.83
4	1.09	0.03	34.5	9.83E-08	454	38.2	9.83E-08	454	0.91	2.76
5	3.25	0.06	26.3	9.71E-08	272	55.2	9.71E-08	272	0.51	0.53
6	11.0	0.07	44.7	9.42E-08	263	84.0	9.44E-08	263	0.59	1.02
7	5.95	0.09	8.84	8.53E-08	43	17.0	8.53E-08	43	0.64	1.68
8	2.94	0.04	0.09	2.10E-08	1.5	1.75	2.10E-08	1.5	0.64	3.54
9	0.11	0.01	53.2	4.24E-03	2000	24.4	4.92E-03	2000	2.18	1.31
10	0.12	0.01	0.16	4.89E-11	4	0.08	6.45E-11	4	1.36	1.18
11	2.89	0.09	0.44	6.10E-09	2.5	0.48	6.10E-09	2.5	1.00	33.2
12	0.36	0.03	36.6	7.05E-08	278.5	10.6	8.82E-08	270.5	3.35	8.04

表5. csrilut(20,10^-7) 预处理CG与BiCGStab方法

		ilut(20,10-7)		CPU		GPU		加速比
#	实际时间(秒)	复制时间(秒)	求解时间(秒)	\(\frac{\left\| \mathbf{r}_{i} \right\|_{2}}{\left\| \mathbf{r}_{0} \right\|_{2}}\)	# 迭代次数	求解时间(秒)	\(\frac{\left\| \mathbf{r}_{i} \right\|_{2}}{\left\| \mathbf{r}_{0} \right\|_{2}}\)	# 迭代次数	对比 ilut (20,10-7)	对比 ilu0
1	0.82	0.02	47.6	9.90E-08	1297	159	9.86E-08	1292	0.30	25.2
2	9.21	0.11	32.1	8.69E-08	193	84.6	8.67E-08	193	0.44	1.16
3	10.04	0.04	6.26	9.64E-08	90	4.75	9.64E-08	90	1.10	2.36
4	8.12	0.10	15.7	9.02E-08	148	22.5	9.02E-08	148	0.78	1.84
5	8.60	0.10	21.2	9.52E-08	158	53.6	9.52E-08	158	0.48	0.54
6	35.2	0.11	29.2	9.88E-08	162	80.5	9.88E-08	162	0.56	1.18
7	23.1	0.14	3.79	7.50E-08	14	12.1	7.50E-08	14	0.76	3.06
8	5.23	0.05	0.14	1.19E-09	1.5	2.37	1.19E-09	1.5	0.70	6.28
9	0.12	0.01	55.1	3.91E-03	2000	24.4	2.27E-03	2000	2.25	1.36
10	0.14	0.01	0.14	9.35E-08	3.5	0.07	7.19E-08	3.5	1.28	1.18
11	218	0.12	0.43	9.80E-08	2	0.66	9.80E-08	2	1.00
12	15.0	0.21	12.2	3.45E-08	31	4.95	3.45E-08	31	1.35	5.93

在GPU上使用不同不完全分解预处理的BiCGStab和CG迭代方法性能总结如图3所示，其中"*"表示该方法未能收敛到所需的容差。需要注意的是，在我们的数值实验中，随着阈值参数放宽和分解变得更加密集，不完全分解的性能通常会下降，这由于稀疏三角求解中行间数据依赖而抑制了并行性。因此，在GPU上获得最佳性能的是零填充的不完全-LU和Cholesky分解，这将作为我们的参考基准。

使用不完全LU Cholesky预处理的BiCGStab与CG算法性能对比

尽管采用更宽松阈值的不完全分解通常更接近精确分解，从而减少迭代步骤，但其计算成本也显著更高。此外需注意，即使迭代步骤数量减少，每一步的计算开销却更大。由于这些权衡因素，在我们的数值实验中，迭代方法的总时间（分解阶段耗时与迭代求解阶段耗时之和）并不一定会随着阈值放宽而减少。

基于GPU上使用csrilu0预处理器和CPU上使用全部四种预处理器的预处理迭代方法总耗时的加速比展示在图4中。值得注意的是，在我们的数值实验中，对于大多数矩阵而言，使用cuSPARSE和cuBLAS库实现的迭代方法确实优于MKL。

GPU上使用csrilu0的预条件BiCGStab和CG加速比 vs CPU上使用全部

最后，获得的加速比平均值显示在图5中，其中我们排除了ilut(10,10^-5)对cage13矩阵的运行结果，以及ilut(20,10^-7)对offshore和cage13矩阵的不完全分解运行结果，因为它们的加速比不成比例。不过，包含这些运行的加速比已在同一图中以括号形式显示。因此我们可以得出结论：在GPU上实现的不完全-LU和Cholesky预处理的BiCGStab及CG迭代方法，相比其CPU实现平均获得了超过2倍的加速。

BiCGStab和CG在GPU（使用csrilu0）与CPU（使用全部）上的平均加速比

4. 总结

迭代方法的性能高度依赖于当前系数矩阵的稀疏模式。在我们的数值实验中，使用cuSPARSE和cuBLAS库在GPU上实现的不完全LU分解和Cholesky预处理的BiCGStab及CG迭代方法，相比其MKL实现平均获得了2倍的加速。

稀疏矩阵向量乘法和三角求解（分为只需执行一次的较慢分析阶段和可多次执行的较快求解阶段）是这些迭代方法的核心构建模块。实际上，所获得的加速效果通常主要受算法求解阶段耗时的影响。

最后需要指出的是，使用多右端项可以增加可用并行度，并能在预处理迭代方法中带来显著的相对性能提升。此外，采用CUDA并行编程范式开发不完全LU分解和Cholesky分解，还能进一步提升所获得的加速效果。

5. 致谢

本白皮书由Maxim Naumov为NVIDIA公司撰写。

允许出于任何用途对本作品的全部或部分内容进行数字或纸质复制，前提是复制件保留本声明及首页的完整引用信息，且无需支付费用。

6. 参考文献

[1] E. Anderson 和 Y. Saad 在并行计算机上求解稀疏三角线性系统，《国际高速计算杂志》，第73-95页，1989年。

[2] R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. van der Vorst，《线性系统求解模板：迭代方法的构建模块》，SIAM出版社，费城，宾夕法尼亚州，1994年。

[3] N. Bell 和 M. Garland，《在面向吞吐量的处理器上实现稀疏矩阵-向量乘法》，高性能计算网络、存储与分析会议论文集(SC09)，ACM出版社，第1-11页，2009年。

[4] A. Greenbaum, 在NYU Ultracomputer原型机上使用带并行扩展的Fortran求解稀疏三角线性方程组, 技术报告99, NYU Ultracomputer笔记, 纽约大学, 纽约州, 1986年4月.

[5] D. B. Kirk 和 W. W. Hwu，《大规模并行处理器编程：实践方法》，Elsevier出版社，2010年。

[6] J. Mayer，《稀疏三角矩阵线性系统求解的并行算法》，计算，第291-312页（86），2009年。

[7] R. Mirchandaney, J. H. Saltz 和 D. Baxter，《运行时循环的并行化与调度》，IEEE计算机汇刊，第(40)期，1991年。

[8] M. Naumov, 预处理迭代方法中稀疏三角线性系统在GPU上的并行求解, NVIDIA技术报告, NVR-2011-001, 2011.

[9] J. Nickolls, I. Buck, M. Garland 和 K. Skadron, 使用CUDA进行可扩展并行编程, Queue, 第40-53页 (6-2), 2008年。

[10] E. Rothberg 和 A. Gupta, 层次存储多处理器上的并行ICCG - 解决三角求解瓶颈问题, Parallel Comput., 第719-741页(18), 1992年.

[11] Y. Saad, 《稀疏线性系统的迭代方法》, SIAM出版社, 费城, 宾夕法尼亚州, 第二版, 2003年。

[12] J. H. Saltz，《多处理器上求解稀疏三角系统的聚合方法》，SIAM科学与统计计算期刊，第123-144页(11)，1990年。

[13] J. Sanders 和 E. Kandrot，《CUDA实战：通用GPU编程入门》，Addison-Wesley出版社，2010年。

[14] M. Wolf, M. Heroux 和 E. Boman，《多线程稀疏三角求解的性能影响因素》，第九届高性能计算与科学国际会议(VECPAR)，2010年。

[15] S. Williams, L. Oliker, R. Vuduc, J. Shalf, K. Yelick 和 J. Demmel，《新兴多核平台上稀疏矩阵向量乘法的优化》，Parallel Comput.，第178-194页（35-3），2009年。

[16] NVIDIA cuSPARSE and cuBLAS Libraries,http://www.nvidia.com/object/cuda_develop.html

[17] 英特尔数学核心函数库，http://software.intel.com/en-us/articles/intel-mkl

[18] 佛罗里达大学稀疏矩阵集合，http://www.cise.ufl.edu/research/sparse/matrices/。

7. 通知

7.1. 注意事项

本文档仅供信息参考之用，不应视为对产品功能、状态或质量的保证。NVIDIA公司（“NVIDIA”）对本文件所含信息的准确性或完整性不作任何明示或暗示的陈述或保证，并对其中可能存在的错误不承担任何责任。NVIDIA对于因使用此类信息而产生的后果、或因使用该信息导致的第三方专利或其他权利侵权概不负责。本文件不构成对开发、发布或交付任何材料（定义见下文）、代码或功能的承诺。

NVIDIA保留随时对本文件进行更正、修改、增强、改进以及任何其他变更的权利，恕不另行通知。

客户在下单前应获取最新的相关信息，并确认这些信息是最新且完整的。

除非NVIDIA与客户授权代表签署的单独销售协议中另有约定，否则NVIDIA产品的销售均以订单确认时提供的NVIDIA标准销售条款和条件为准（以下简称"销售条款"）。NVIDIA特此明确反对将任何客户通用条款适用于本文件所述NVIDIA产品的采购。本文件不直接或间接构成任何合同义务。

NVIDIA产品并非设计、授权或保证适用于医疗、军事、航空、航天或生命支持设备，也不适用于那些可以合理预期NVIDIA产品故障或失灵会导致人身伤害、死亡、财产或环境损害的应用场景。NVIDIA对于在此类设备或应用中使用和/或包含NVIDIA产品不承担任何责任，因此客户需自行承担相关风险。

NVIDIA不声明或保证基于本文档的产品适用于任何特定用途。NVIDIA未必会对每个产品的所有参数进行测试。客户应全权负责评估和确定本文档所含信息的适用性，确保产品适合并满足客户计划的应用需求，并执行必要的应用测试以避免应用或产品出现故障。客户产品设计中的缺陷可能会影响NVIDIA产品的质量和可靠性，并可能导致超出本文档范围的其他或不同的条件和/或要求。对于任何因以下原因导致的故障、损坏、成本或问题，NVIDIA不承担任何责任：(i) 以违反本文档的任何方式使用NVIDIA产品或(ii) 客户产品设计。

本文档不授予任何NVIDIA专利权、版权或其他NVIDIA知识产权的明示或暗示许可。NVIDIA发布的关于第三方产品或服务的信息，不构成NVIDIA对这些产品或服务的使用许可或担保认可。使用此类信息可能需要获得第三方基于其专利或其他知识产权的许可，或需要获得NVIDIA基于其专利或其他知识产权的许可。

本文件中的信息仅可在获得NVIDIA事先书面批准、未经改动完整复制且完全符合所有适用的出口法律法规，并附带所有相关条件、限制和声明的情况下进行复制。

本文件及所有NVIDIA设计规格、参考板、文件、图纸、诊断工具、清单和其他文档（统称及单独称为"材料"）均以"现状"提供。NVIDIA不对材料作出任何明示或默示的保证，包括但不限于对不侵权、适销性和特定用途适用性的默示保证免责。在法律允许的最大范围内，NVIDIA不就因使用本文件导致的任何损害承担责任，包括但不限于任何直接、间接、特殊、附带、惩罚性或后果性损害，无论损害成因如何，也无论责任理论为何，即使NVIDIA已被告知发生此类损害的可能性。不论客户因任何原因可能遭受的任何损害，NVIDIA对客户就本文所述产品的全部及累计责任应受产品销售条款的限制。

7.2. OpenCL

OpenCL是苹果公司的商标，经Khronos Group Inc.授权使用。

7.3. 商标

NVIDIA和NVIDIA标识是美国及其他国家NVIDIA公司的商标或注册商标。其他公司及产品名称可能是其各自关联公司的商标。