我有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传入，并被解决方案覆盖。 Lp，Li和Lx分别是稀疏矩阵的标准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/benchmark和this 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/