使用 OpenACC 正确实现稀疏矩阵的三角求解器

Question

我目前正在为稀疏矩阵实现三角求解器，并尝试使用 OpenACC 指令进行加速。鉴于我的矩阵因子 LU 采用稀疏 CSR 格式，OpenACC 已经成功地解决了 L 因子，但与应用程序的真实解决方案相比，U 因子给出了完全错误的结果。下面是用于反向替换任务的加速内核代码：

#pragma acc kernels deviceptr( ia, ja, factorValsU, y, x )
{
    for ( int idx = size; idx > 0; idx-- )
    {
        double temp = 0.0;
        int rowInit = ia[ idx - 1];
        int rowEnd  = ia[ idx ];
        #pragma acc loop vector reduction( + : temp)
        for ( int k = rowInit + 1; k < rowEnd; k++ )
        {
            temp += factorValsU[ k ] * x[ ja[ k ] ];
        }
        x[ idx ] = (y[ idx ] - temp) / factorValsU[ rowInit ];
    }
 }

我不知道为什么这个内核会产生不正确的结果。我已经尝试了一个不同版本的内核，其中矩阵向后保存，即从下到上，原则上可以通过以下内核解决：

#pragma acc kernels deviceptr( ia, ja, factorValsU, y, x )
{
    for ( int idx = 0; idx < size; idx++ )
    {
        double temp = 0.0;
        int rowInit = ia[ idx ];
        int rowEnd  = ia[ idx + 1 ];
        #pragma acc loop vector reduction( + : temp)
        for ( int k = rowInit + 1; k < rowEnd; k++ )
        {
            temp += factorValsU[ k ] * x[ ja[ k ] ];
        }
        x[ size - idx ] = (y[ size - idx ] - temp) / factorValsU[ rowInit ];
    }
 }

但结果总是错误的。我是否错过了一些关于使用 OpenACC 指令装饰常规代码以获得正确结果的基本知识？

如前所述，L 因子的前向替换工作正常，为了完整起见，我在此处发布代码。

#pragma acc kernels deviceptr( ia, ja, factorValsL, y, x )
{
    for ( int idx = 0; idx < size; idx++ )
    {
        double temp = 0.0;
        int rowInit = ia[ idx ];
        int rowEnd  = ia[ idx + 1 ];
        #pragma acc loop vector reduction( + : temp)
        for ( int k = rowInit; k < rowEnd; k++ )
        {
            temp += factorValsL[ k ] * x[ ja[ k ] ];
        }
        x[ idx ] = y[ idx ] - temp;
    }
 }

注意前向替换（有效）和后向替换（均无效）的内核之间的细微差别，是保存结果的内存区域：

x[ idx ] = y[ idx ] - temp     for the L factor 
x[ size - idx ] = (y[ size - idx ] - temp) / factorValsU[ rowInit ] for the U factor;

U 因子求解器是否有某些原因会计算错误结果，原因是内存中分配（和讲座）的顺序？

为了完整起见，pgi18.4编译器提供的关于内核的信息是：

triangularSolverU_acc(int, const int *, const int *, const double *, const double *, double *, bool):
614, Complex loop carried dependence of y->,x->,factorVals-> prevents parallelization
     Loop carried dependence of x-> prevents parallelization
     Loop carried backward dependence of x-> prevents vectorization
     Accelerator kernel generated
     Generating Tesla code
    614, #pragma acc loop seq
    621, #pragma acc loop vector(128) /* threadIdx.x */
         Generating reduction(+:temp)
621, Loop is parallelizable

这说明外循环已经序列化，内循环是一个reduction。

Answer 1

对于“内核”，编译器必须证明循环不包含任何依赖项，因此可以安全并行化。但是，由于您的代码包含指针，并且指针可能别名为同一内存，因此编译器无法证明这一点，因此安全的事情也可以使循环按顺序运行。要覆盖编译器分析，您可以在外部 for 循环之前添加“#pragma acc loopdependent”。 “独立”是对编译器的断言，即循环可以安全并行化。或者，您可以使用“并行循环”指令代替“内核”，因为“并行”意味着“独立”。

对于错误的答案，这通常是由于主机和设备副本之间的数据未正确同步。您如何管理数据移动？由于您使用的是“deviceptr”，这意味着您使用的是 CUDA。

另外，如果您可以发布一个完整的复制示例，将更容易帮助确定问题。