c++ - 单/ double SpMV 在 CPU 上的性能

Question

由于非常低的算术强度，稀疏矩阵 vector 乘积是一种内存限制操作。由于浮点存储格式需要 4+4=8 个字节的每个非零，而 double (值和列索引)需要 4+8=12 个字节，因此当切换到浮点时，执行速度应该可以提高 33%。我构建了一个基准，它组装了一个 1000000x1000000 矩阵，每行有 200 个非零值，然后从 20 次乘法中取最小值。 github上的源代码here .

结果和我预想的差不多。当我在我的英特尔酷睿 i7-2620M 上运行基准测试时，我发现执行速度提高了 30%。在 spec 中的 21.3 GB/s 中，可以看出带宽从大约 19.0 GB/s( double )下降到大约 18.0 GB/s(浮点)的微小差异。 .

现在，由于矩阵的数据比 vector 的数据大将近 1000，因此对于只有矩阵是单精度但 vector 保持为 double 。我继续尝试这个，然后确保使用较低的精度进行计算。然而，当我运行它时，有效带宽使用率突然下降到大约 14.4 GB/s，执行速度只比完整的双版本快 12%。如何理解这一点？

---

我使用的是 Ubuntu 14.04 和 GCC 4.9.3。

Run times:

// double(mat)-double(vec)
Wall time: 0.127577 s
Bandwidth: 18.968 GB/s
Compute:   3.12736 Gflop/s

// float(mat)-float(vec)
Wall time: 0.089386 s
Bandwidth: 18.0333 GB/s
Compute:   4.46356 Gflop/s

// float(mat)-double(vec)
Wall time: 0.112134 s
Bandwidth: 14.4463 GB/s
Compute:   3.55807 Gflop/s

---

更新

请参阅下面 Peter Cordes 的回答。简而言之，从 double 到浮点转换的循环迭代之间的依赖关系是造成开销的原因。通过展开循环(参见 github 上的 unroll-loop 分支)，float-double 和 float-float 版本都重新获得了全部带宽使用!

New run times:

// float(mat)-float(vec)
Wall time: 0.084455 s
Bandwidth: 19.0861 GB/s
Compute:   4.72417 Gflop/s

// float(mat)-double(vec)
Wall time: 0.0865598 s
Bandwidth: 18.7145 GB/s
Compute:   4.6093 Gflop/s

Answer 1

必须即时转换的双浮点循环无法发出那么快。随着一些循环展开，gcc 可能会做得更好。

你的 i7-2620M is a dual-core with hyperthreading .当瓶颈是 CPU uop 吞吐量，而不是分支预测错误、缓存未命中，甚至只是长延迟链时，超线程无济于事。仅通过标量运算使内存带宽饱和并不容易。

---

来自 Godbolt Compiler Explorer 上代码的 asm 输出: gcc 5.3 做了大致相同的内部循环，顺便说一句，所以在这种情况下，使用旧的 gcc 版本并不会损失太多。

double-double version inner loop (gcc 4.9.3 -O3 -march=sandybridge -fopenmp):

## inner loop of <double,double>mult() with fused-domain uop counts
.L7:
    mov     edx, eax  # 1 uop
    add     eax, 1    # 1 uop
    mov     ecx, DWORD PTR [r9+rdx*4] # 1 uop
    vmovsd  xmm0, QWORD PTR [r10+rdx*8] # 1 uop
    vmulsd  xmm0, xmm0, QWORD PTR [r8+rcx*8]    # 2 uops
    vaddsd  xmm1, xmm1, xmm0    # 1 uop
    cmp     eax, esi  #  (macro-fused)
    jne     .L7       #  1 uop

总计:8 个融合域微指令，每两个时钟可以发出一个迭代器。它还可以执行得那么快:其中三个微指令是负载，SnB 每 2 个时钟可以执行 4 次负载。剩下 5 个 ALU 微指令(因为 SnB 无法在重命名阶段消除 reg-reg 移动，这是在 IvB 中引入的)。反正单个执行端口没有明显的瓶颈。 SnB 的三个 ALU 端口每两个周期最多可以处理六个 ALU 微指令。

没有微融合 because of using two-register addressing modes .

---

双浮点型内循环:

## inner loop of <double,float>mult() with fused-domain uop counts
.L7:
    mov     edx, eax  # 1 uop
    vxorpd  xmm0, xmm0, xmm0    # 1 uop (no execution unit needed).
    add     eax, 1    # 1 uop
    vcvtss2sd       xmm0, xmm0, DWORD PTR [r9+rdx*4]      # 2 uops
    mov     edx, DWORD PTR [r8+rdx*4] # 1 uop
    vmulsd  xmm0, xmm0, QWORD PTR [rsi+rdx*8]   # 2 uops
    vaddsd  xmm1, xmm1, xmm0    # 1 uop
    cmp     eax, ecx  # (macro-fused)
    jne     .L7       # 1 uop

gcc 使用 xorpd 来打破循环携带的依赖链。 cvtss2sd 错误地依赖于 xmm0 的旧值，因为它的设计很糟糕，并且没有将寄存器的上半部分归零。 (movsd 当用作加载时为零，但当用作 reg-reg 移动时则不会。在这种情况下，请使用 movaps 除非你想合并。)

因此，10 个融合域微指令:可以每三个时钟迭代一次。我认为这是唯一的瓶颈，因为它只是一个需要执行端口的额外 ALU uop。 (SnB handles zeroing-idioms in the rename stage, so xorpd doesn't need one)。 cvtss2sd 是一条 2 uop 指令，即使 gcc 使用单寄存器寻址模式，它显然也无法微融合。它的吞吐量为每个时钟一个。 (在 Haswell 上，当与寄存器 src 和 dest 一起使用时，它是一个 2 uop 指令，而在 Skylake 上，吞吐量减少到每 2 个时钟一个，根据 Agner Fog's testing 。)这仍然不是这个循环的瓶颈不过，天湖。 Haswell/Skylake 上仍然有 10 个融合域 uops，这仍然是瓶颈。

---

-funroll-loops 应该有助于 gcc 4.9.3

gcc 做得还不错，代码如下

    mov     edx, DWORD PTR [rsi+r14*4]        # D.56355, *_40
    lea     r14d, [rax+2]     # D.56355,
    vcvtss2sd       xmm6, xmm4, DWORD PTR [r8+r14*4]      # D.56358, D.56358, *_36
    vmulsd  xmm2, xmm1, QWORD PTR [rcx+rdx*8]   # D.56358, D.56358, *_45
    vaddsd  xmm14, xmm0, xmm13  # tmp, tmp, D.56358
    vxorpd  xmm1, xmm1, xmm1    # D.56358

    mov     edx, DWORD PTR [rsi+r14*4]        # D.56355, *_40
    lea     r14d, [rax+3]     # D.56355,
    vcvtss2sd       xmm10, xmm9, DWORD PTR [r8+r14*4]     # D.56358, D.56358, *_36
    vmulsd  xmm7, xmm6, QWORD PTR [rcx+rdx*8]   # D.56358, D.56358, *_45
    vaddsd  xmm3, xmm14, xmm2   # tmp, tmp, D.56358
    vxorpd  xmm6, xmm6, xmm6    # D.56358

如果没有循环开销，每个元素的工作减少到 8 个融合域微指令，并且它不是一个每 3 个周期只发出 2 微指令的小循环(因为 10 不是 4 的倍数)。

它可以通过使用位移来保存 lea 指令，例如[r8+rax*4 + 12]。 IDK 为什么 gcc 选择不这样做。

甚至 -ffast-math 都无法对其进行矢量化。这可能没有意义，因为从稀疏矩阵进行收集会超过从非稀疏 vector 加载 4 或 8 个连续值的好处。 (内存中的 insertps 是一个 2-uop 指令，即使在单寄存器寻址模式下也不能微融合。)

在 Broadwell 或 Skylake 上，vgatherdps 可能足够快以提供加速。可能是 Skylake 的一大加速。 (可以收集 8 个单精度 float ，吞吐量为每 5 个时钟 8 个 float 。或者 vgatherqpd 可以收集 4 个 double float ，吞吐量为每 4 个时钟 4 个 double float )。这将为您设置一个 256b vector FMA。