我有n x n个下三角矩阵A的压缩稀疏列(csc)表示,在主对角线上有零,并想解​​决in中的b

(A + I)' * x = b

这是我计算此例程的例程:
void backsolve(const int*__restrict__ Lp,
               const int*__restrict__ Li,
               const double*__restrict__ Lx,
               const int n,
               double*__restrict__ x) {
  for (int i=n-1; i>=0; --i) {
      for (int j=Lp[i]; j<Lp[i+1]; ++j) {
          x[i] -= Lx[j] * x[Li[j]];
      }
  }
}

因此,b通过参数x传入,并被解决方案覆盖。 LpLiLx分别是稀疏矩阵的标准csc​​表示形式中的行,索引和数据指针。该功能是程序中 HitTest 门的地方,与
x[i] -= Lx[j] * x[Li[j]];

是花费的大部分时间。用gcc-8.3 -O3 -mfma -mavx -mavx512f编译给出
backsolve(int const*, int const*, double const*, int, double*):
        lea     eax, [rcx-1]
        movsx   r11, eax
        lea     r9, [r8+r11*8]
        test    eax, eax
        js      .L9
.L5:
        movsx   rax, DWORD PTR [rdi+r11*4]
        mov     r10d, DWORD PTR [rdi+4+r11*4]
        cmp     eax, r10d
        jge     .L6
        vmovsd  xmm0, QWORD PTR [r9]
.L7:
        movsx   rcx, DWORD PTR [rsi+rax*4]
        vmovsd  xmm1, QWORD PTR [rdx+rax*8]
        add     rax, 1
        vfnmadd231sd    xmm0, xmm1, QWORD PTR [r8+rcx*8]
        vmovsd  QWORD PTR [r9], xmm0
        cmp     r10d, eax
        jg      .L7
.L6:
        sub     r11, 1
        sub     r9, 8
        test    r11d, r11d
        jns     .L5
        ret
.L9:
        ret

根据vtune,
vmovsd  QWORD PTR [r9], xmm0

是最慢的部分。我几乎没有组装方面的经验,并且对如何进一步诊断或优化此操作一无所知。我尝试使用不同的标志进行编译以启用/禁用SSE,FMA等,但是没有任何效果。

处理器:Xeon Skylake

问题我应该怎么做才能优化此功能?

最佳答案

这应该在很大程度上取决于矩阵和所使用平台的精确稀疏性模式。我在大小为3006的this matrix的下三角上使用gcc 8.3.0和编译器标志-O3 -march=native(在我的CPU上为-march=skylake)测试了几件事,带有19554个非零条目。希望这在某种程度上接近您的设置,但是无论如何我希望它们能使您对从哪里开始有所了解。

为了计时,我使用了google/benchmarkthis source file。它定义了对问题中给定的实现进行基准测试的benchBacksolveBaseline和对所提出的“优化”实现进行了基准测试的benchBacksolveOptimized。还有benchFillRhs分别对这两个函数使用的基准进行基准测试,以为右侧生成一些不完全平凡的值。为了获得“纯”反解析的时间,应减去benchFillRhs花费的时间。

1.严格向后迭代

实现中的外循环向后遍历各列,而内循环向前遍历当前各列。似乎也可以向后遍历每一列更加一致:

for (int i=n-1; i>=0; --i) {
    for (int j=Lp[i+1]-1; j>=Lp[i]; --j) {
        x[i] -= Lx[j] * x[Li[j]];
    }
}

这几乎不会更改程序集(https://godbolt.org/z/CBZAT5),但是基准测试时间显示出可衡量的改进:
------------------------------------------------------------------
Benchmark                        Time             CPU   Iterations
------------------------------------------------------------------
benchFillRhs                  2737 ns         2734 ns      5120000
benchBacksolveBaseline       17412 ns        17421 ns       829630
benchBacksolveOptimized      16046 ns        16040 ns       853333

我认为这是由于某种程度上更可预测的缓存访问引起的,但我没有对其进行深入研究。

2.内循环中的负载/存储量更少

由于A是较低的三角形,因此我们有i < Li[j]。因此,我们知道x[Li[j]]不会由于内部循环中x[i]的更改而更改。我们可以使用一个临时变量将这些知识放入我们的实现中:
for (int i=n-1; i>=0; --i) {
    double xi_temp = x[i];
    for (int j=Lp[i+1]-1; j>=Lp[i]; --j) {
        xi_temp -= Lx[j] * x[Li[j]];
    }
    x[i] = xi_temp;
}

这使得gcc 8.3.0将存储从内部循环内部移动到内存中,直到其结束(https://godbolt.org/z/vM4gPD)之后。我的系统上测试矩阵的基准显示出一点改进:
------------------------------------------------------------------
Benchmark                        Time             CPU   Iterations
------------------------------------------------------------------
benchFillRhs                  2737 ns         2740 ns      5120000
benchBacksolveBaseline       17410 ns        17418 ns       814545
benchBacksolveOptimized      15155 ns        15147 ns       887129

3.展开循环

尽管在首次建议的代码更改之后clang已经开始展开循环,但是gcc 8.3.0仍然没有。因此,让我们尝试通过额外传递-funroll-loops
------------------------------------------------------------------
Benchmark                        Time             CPU   Iterations
------------------------------------------------------------------
benchFillRhs                  2733 ns         2734 ns      5120000
benchBacksolveBaseline       15079 ns        15081 ns       953191
benchBacksolveOptimized      14392 ns        14385 ns       963441

请注意,基线也有所改善,因为该实现中的循环也已展开。我们的优化版本还可以从循环展开中受益,但可能没有我们喜欢的那么多。查看生成的程序集(https://godbolt.org/z/_LJC5f),似乎gcc在展开8次时可能有点过头了。对于我的设置,实际上我可以通过一个简单的手动展开来做得更好。因此,再次删除标记-funroll-loops,并使用以下方式实现展开:
for (int i=n-1; i>=0; --i) {
    const int col_begin = Lp[i];
    const int col_end = Lp[i+1];
    const bool is_col_nnz_odd = (col_end - col_begin) & 1;
    double xi_temp = x[i];
    int j = col_end - 1;
    if (is_col_nnz_odd) {
        xi_temp -= Lx[j] * x[Li[j]];
        --j;
    }
    for (; j >= col_begin; j -= 2) {
        xi_temp -= Lx[j - 0] * x[Li[j - 0]] +
                   Lx[j - 1] * x[Li[j - 1]];
    }
    x[i] = xi_temp;
}

这样,我可以测量:
------------------------------------------------------------------
Benchmark                        Time             CPU   Iterations
------------------------------------------------------------------
benchFillRhs                  2728 ns         2729 ns      5090909
benchBacksolveBaseline       17451 ns        17449 ns       822018
benchBacksolveOptimized      13440 ns        13443 ns      1018182

其他算法

所有这些版本仍然对稀疏矩阵结构使用相同的向后求解的简单实现。从本质上讲,在像这样的稀疏矩阵结构上进行操作可能会给内存流量带来重大问题。至少对于矩阵分解而言,存在更为复杂的方法,这些方法可用于由稀疏结构组装而成的密集子矩阵。示例是超结点方法和多面方法。我对此感到有点模糊,但是我认为这样的方法还将这种思想应用于布局,并使用密集矩阵运算进行下三角向后求解(例如,对于Cholesky型分解)。因此,如果您不被迫坚持直接在稀疏结构上运行的简单方法,那么可能值得研究这类方法。例如,请参阅戴维斯(Davis)的this survey

关于c++ - 优化稀疏下三角线性系统的后向求解,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/60232977/

10-13 03:26