cusparse csrsv_analysis 的性能非常慢

Question

我编写了一个带有 LU 预处理的共轭梯度求解器（用于线性方程组），我使用了 Maxim Naumov 博士在 nvidia 研究社区上的 papers 作为指导，残差更新步骤，这需要求解下三角矩阵系统然后求解一个上三角矩阵系统分为两个阶段：

1. 分析阶段（利用稀疏模式并决定并行化级别）。
2. 解决阶段本身。

根据与该主题相关的所有帖子以及 Naumov 的论文本身，分析阶段比求解阶段慢得多，但它只执行一次，因此在考虑整个执行时间时应该不是问题，但是，在我的程序中，分析阶段需要大约 35-45% 的整个求解时间（执行所有迭代所需的时间！），这很烦人，另一件事是我很确定矩阵的稀疏模式不允许大量并行化，因为几乎所有元素都需要知道先前的元素（因为在 CFD 中，每个节点都需要至少相邻的 6 个节点（六面体体积）围绕它，并且每一行都重复），所以无论如何，分析阶段不会很有用！

ma​​trixLU 这段代码中同时包含上三角矩阵和下三角矩阵，上三角矩阵使用原矩阵对角线，下三角矩阵有单位对角线（LU分解）。

// z = inv(matrixLU)*r
cusparseMatDescr_t descrL = 0 ;
cusparseMatDescr_t descrU = 0 ;
cusparseStatus = cusparseCreateMatDescr(&descrL) ;
cusparseStatus = cusparseCreateMatDescr(&descrU) ;

cusparseSetMatType(descrL,CUSPARSE_MATRIX_TYPE_GENERAL) ;
cusparseSetMatIndexBase(descrL,CUSPARSE_INDEX_BASE_ONE) ;
cusparseSetMatDiagType(descrL,CUSPARSE_DIAG_TYPE_UNIT) ;
cusparseSetMatFillMode(descrL,CUSPARSE_FILL_MODE_LOWER) ;

cusparseSetMatType(descrU,CUSPARSE_MATRIX_TYPE_GENERAL) ;
cusparseSetMatIndexBase(descrU,CUSPARSE_INDEX_BASE_ONE) ;
cusparseSetMatDiagType(descrU,CUSPARSE_DIAG_TYPE_NON_UNIT) ;
cusparseSetMatFillMode(descrU,CUSPARSE_FILL_MODE_UPPER) ;

cusparseSolveAnalysisInfo_t inforL = 0 ;
cusparseSolveAnalysisInfo_t inforU = 0 ;
cusparseStatus = cusparseCreateSolveAnalysisInfo(&inforL) ;
cusparseStatus = cusparseCreateSolveAnalysisInfo(&inforU) ;

double startFA = omp_get_wtime() ;
cusparseStatus = cusparseDcsrsv_analysis(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, NZ, descrL, matrixLU, iRow, jCol, inforL) ;
if(cusparseStatus != CUSPARSE_STATUS_SUCCESS) printf("%s \n\n","cusparseDcsrsv_analysis1 Error !") ;
cusparseStatus = cusparseDcsrsv_analysis(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, NZ, descrU, matrixLU, iRow, jCol, inforU) ;
if(cusparseStatus != CUSPARSE_STATUS_SUCCESS) printf("%s \n\n","cusparseDcsrsv_analysis2 Error !") ;
double finishFA = omp_get_wtime() ;

那么，有人知道为什么分析阶段如此缓慢吗？以及如何加速？ （分析阶段执行时间）/（求解阶段执行时间）比率取决于 GPU 吗？

编辑： 我在很多情况下都尝试了这个求解器，结果很接近，但在我关心的具体情况下，它具有以下条件：

• 大小（N）：~860,000 * 860,000
• 非零数 (NZ)：~6,000,000
• 收敛所需的迭代次数：10
• 分析阶段执行时间：210 ms
• 解决阶段执行时间（总结所有迭代）：350 ms
• 所有浮点运算均以双精度格式执行
• GPU：GeForce GTX 550 Ti
• 操作系统：Windows 7 终极版，64 位

Answer 1

我联系了 NVIDIA 的线性代数库团队，他们提供了以下反馈。

1. 关于性能结果：

   1. CG 迭代法只执行 10 次迭代，直到找到线性系统的解。在这种情况下，这似乎不足以分摊分析阶段消耗的时间。 （编辑：） 收敛到解决方案所需的步骤数因应用程序而异。在一些无条件的迭代方法中不收敛（或需要数千次迭代），而有条件的迭代方法需要数十次（或数百次）迭代才能收敛到一个解决方案。在这个特定应用中，迭代方法在 10 次迭代中收敛这一事实不一定代表所有应用。

   2. 分析和求解阶段的比率为 6 = 210/35，这与我们在其他矩阵上的实验一致，通常在区间 [4,10] 内。目前尚不清楚分析和求解阶段的时间与 CPU 上稀疏三角求解所花费的时间相比如何。 （这将是有趣的信息。）

2. 关于算法：
   1. 不幸的是，您只能采取一些措施来恢复一点性能（没有真正的方法可以加快分析阶段）。例如，您可以将 matrixLU 拆分为单独的下三角部分和上三角部分。使用单独的三角形部分进行计算通常比使用嵌入在一般矩阵中的三角形部分进行计算要快。如果您使用填充为 0 的不完全分解作为预处理器，您还可以利用 GPU 上填充为 0 的不完整 LU/Cholesky，该功能最近在 CUSPARSE 库中可用。 （编辑：）CUDA Toolkit 5.0 版本中提供了填充为 0 的不完全 LU 分解。目前此版本的抢先体验版可以通过 NVIDIA 网站上的registered developers 获得。
   2. 我们从未研究过不同 GPU 之间分析和求解阶段的比率如何变化；我们所做的所有实验都是在 C2050 上进行的。
   3. 分析阶段很慢，因为它必须收集有关哪些行可以一起处理成级别的信息。求解阶段更快，因为它只处理关卡（参见this paper）。

最后，您还提到您认为问题的并行度不够。但是仅从矩阵中很难猜出并行度的数量。如果确实很少有并行性（每一行都依赖于前一行），那么 CUSPARSE 稀疏三角求解可能不是解决这个问题的正确算法。此外，您可以尝试在不进行预处理或使用其他类型的可能更适合您的问题的预处理器的情况下运行 CG 迭代方法。

如前所述，了解这段代码在 GPU 和 CPU 上的性能比较会很有趣。

编辑： 关于一些 cmets... 如前所述，预条件迭代方法的最终迭代次数因不同的应用程序而异。在某些情况下，它可能非常小（例如，在您的应用程序中为 10），但在其他情况下，它可能相对较大（以 100 秒为单位）。这篇论文关注的不是“胜利”，而是算法的权衡。它试图让用户对例程有更深入的了解，以便有效地使用它。同样，如果手头的稀疏模式没有任何并行性，那么这不是适合您的算法。

CPU 和 GPU 的对比，有一些问题（比如密集矩阵-矩阵或稀疏矩阵-向量乘法）是 GPU 擅长的，还有一些问题（比如遍历链表或执行完全顺序代码），它不是。稀疏三角线性系统的解决方案介于这些算法之间（取决于系数矩阵的稀疏模式和线性系统的其他属性，它可能适合您，也可能不适合您）。此外，使用 GPU 的决定不仅仅基于单一算法，例如稀疏三角求解，而是基于应用程序使用的整套算法。最终，GPU 在加速和提高整个应用程序的性能方面通常非常成功。

最后，关于 trsm 与 csrsv，我们并不感到惊讶的是，在密集存储中执行操作的小矩阵比在稀疏存储中更快（即使这些小矩阵可能是稀疏的）。这通常不仅适用于三角求解，而且适用于其他线性代数运算（不过，根据运算和矩阵的稀疏模式，它可能发生在不同的交叉点）。同样，稀疏三角求解算法被设计为在迭代设置中运行，其中分析完成一次​​，求解完成多次。因此，运行一次（并在分析阶段计数）会导致更高的交叉点，从而使此操作在稀疏而不是密集格式中更有效，这并不奇怪。