cusparse csrsv_analysis的性能非常慢

Question

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

1. 分析阶段(利用稀疏性模式并决定并行化级别)
2. 溶液相本身

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

matrixLU 在这个代码中包含上下三角矩阵，上三角矩阵使用原始矩阵对角线，下三角矩阵具有单位对角线(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 nn","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 nn","cusparseDcsrsv_analysis2 Error !") ;
double finishFA = omp_get_wtime() ;

所以，有人知道为什么分析阶段如此缓慢吗？以及如何加速？(分析阶段执行时间)/(求解阶段执行时间？

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

• 大小(N)：~860000*860000
• 非零数(NZ)：~6000000
• 收敛所需的迭代次数：10
• 分析阶段执行时间：210毫秒
• 求解阶段执行时间(所有迭代的总和)：350毫秒
• 所有浮点运算都以双精度格式执行
• 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网站上获得此版本的早期访问版本
   2. 我们从未研究过分析和求解阶段的比率如何在不同的GPU中变化；我们所有的实验都是在C2050上进行的
   3. 分析阶段很慢，因为它必须收集关于哪些行可以一起处理为级别的信息。求解阶段更快，因为它只处理级别(参见本文)

最后，您还提到您认为问题中没有足够的并行性。但是，仅仅从矩阵中就很难猜测并行性的大小。如果确实存在很少的并行性(每一行都取决于前一行)，那么CUSPARSE稀疏三角形求解可能不是解决这个问题的正确算法。此外，您可以尝试在没有预处理的情况下运行CG迭代方法，或者使用其他类型的预处理程序，这些预处理程序可能更适合您的问题。

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

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

关于CPU和GPU的比较，有一些问题(例如密集矩阵矩阵或稀疏矩阵向量乘法)GPU很好，还有一些问题(如遍历链表或执行完全顺序的代码)则不然。稀疏三角形线性系统的解决方案介于这些算法之间(取决于系数矩阵的稀疏性模式和线性系统的其他特性，它可能对您有效或无效)。此外，使用GPU的决定不仅仅基于单个算法，例如稀疏三角形求解，而是基于应用程序使用的整个算法集。最终，GPU在加速和提高整个应用程序的性能方面往往非常成功。

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