single/double 精度 SpMV 在 CPU 上的性能
Performance of single/double precision SpMV on CPUs
由于运算强度非常低，稀疏矩阵向量乘积是一种内存受限操作。由于浮点存储格式需要 4+4=8 个字节的每个非零，而双精度数（值和列索引）需要 4+8=12 个字节，因此当切换到浮点时，执行速度应该可以提高 33%。我构建了一个基准，它组装了一个 1000000x1000000 矩阵，每行有 200 个非零值，然后从 20 次乘法中取最小值。 github here.
上的源代码
结果和我预想的差不多。当我 运行 在我的英特尔酷睿 i7-2620M 上进行基准测试时，我发现执行速度提高了 30%。在 spec 中的 21.3 GB/s 中，可以看出带宽从大约 19.0 GB/s（双倍）下降到大约 18.0 GB/s（浮点）的微小差异。 =15=]
现在，由于矩阵的数据比向量的数据大将近 1000，因此对于只有矩阵是单精度的情况，人们会期望获得更快的性能，但是向量仍然是双打。我继续尝试这个，然后确保使用较低的精度进行计算。然而，当我 运行 它时，有效带宽使用率突然下降到大约 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 的回答。简而言之，从双精度到浮点转换的循环迭代之间的依赖关系是造成开销的原因。通过展开循环（参见 github 处的展开循环分支），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
必须即时转换的双浮点循环不能发出那么快。随着一些循环展开，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。 (). cvtss2sd is a 2 uop instruction that apparently can't micro-fuse even if gcc used a one-register addressing mode. It has a throughput of one per clock. (On Haswell, it's a 2 uop instruction when used with a register src and dest, and on Skylake the throughput is reduced to one per 2 clocks, according to Agner Fog's testing.) 不过，这仍然不是 Skylake 上此循环的瓶颈。 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 都无法对其进行矢量化。这可能没有意义，因为从稀疏矩阵进行收集会超过从非稀疏向量加载 4 或 8 个连续值的好处。 （内存中的 insertps 是一条 2-uop 指令，即使使用单寄存器寻址模式也无法微融合。）
在 Broadwell 或 Skylake 上，vgatherdps 可能足够快以提供加速。可能是 Skylake 的一大加速。 （可以收集 8 个单精度浮点数，吞吐量为每 5 个时钟 8 个浮点数。或者 vgatherqpd 可以收集 4 个双精度浮点数，吞吐量为每 4 个时钟 4 个双精度浮点数）。这将为您设置 256b 矢量 FMA。
由于运算强度非常低，稀疏矩阵向量乘积是一种内存受限操作。由于浮点存储格式需要 4+4=8 个字节的每个非零，而双精度数（值和列索引）需要 4+8=12 个字节，因此当切换到浮点时，执行速度应该可以提高 33%。我构建了一个基准，它组装了一个 1000000x1000000 矩阵，每行有 200 个非零值，然后从 20 次乘法中取最小值。 github here.
上的源代码结果和我预想的差不多。当我 运行 在我的英特尔酷睿 i7-2620M 上进行基准测试时，我发现执行速度提高了 30%。在 spec 中的 21.3 GB/s 中，可以看出带宽从大约 19.0 GB/s（双倍）下降到大约 18.0 GB/s（浮点）的微小差异。 =15=]
现在，由于矩阵的数据比向量的数据大将近 1000，因此对于只有矩阵是单精度的情况，人们会期望获得更快的性能，但是向量仍然是双打。我继续尝试这个，然后确保使用较低的精度进行计算。然而，当我 运行 它时，有效带宽使用率突然下降到大约 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 的回答。简而言之，从双精度到浮点转换的循环迭代之间的依赖关系是造成开销的原因。通过展开循环（参见 github 处的展开循环分支），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
必须即时转换的双浮点循环不能发出那么快。随着一些循环展开，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。 (cvtss2sd is a 2 uop instruction that apparently can't micro-fuse even if gcc used a one-register addressing mode. It has a throughput of one per clock. (On Haswell, it's a 2 uop instruction when used with a register src and dest, and on Skylake the throughput is reduced to one per 2 clocks, according to Agner Fog's testing.) 不过，这仍然不是 Skylake 上此循环的瓶颈。 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 都无法对其进行矢量化。这可能没有意义，因为从稀疏矩阵进行收集会超过从非稀疏向量加载 4 或 8 个连续值的好处。 （内存中的 insertps 是一条 2-uop 指令，即使使用单寄存器寻址模式也无法微融合。）
在 Broadwell 或 Skylake 上，vgatherdps 可能足够快以提供加速。可能是 Skylake 的一大加速。 （可以收集 8 个单精度浮点数，吞吐量为每 5 个时钟 8 个浮点数。或者 vgatherqpd 可以收集 4 个双精度浮点数，吞吐量为每 4 个时钟 4 个双精度浮点数）。这将为您设置 256b 矢量 FMA。