稀疏矩阵向量乘积由于非常低的算术强度而成为内存绑定操作。由于浮动存储格式将需要每个非零4 + 4 = 8个字节,而双打(值和列索引)则需要4 + 8 = 12个字节,因此切换到浮动类型时,应该可以预期将执行速度提高约33%。我构建了一个基准,该基准将一个1000000x1000000矩阵组装成每行200个非零,然后从20个乘法中取最小值。 github here上的源代码。
结果大致符合我的预期。当我在Intel Core i7-2620M上运行基准测试时,执行速度提高了30%。在spec中的21.3 GB / s中,带宽下降从大约19.0 GB / s(两倍)下降到大约18.0 GB / s(浮动),可以看到很小的差异。
现在,由于矩阵的数据比矢量的数据大将近1000,因此人们希望对于仅矩阵为单精度但矢量保持为双精度的情况也应获得更快的性能。我继续尝试,然后确保使用较低的精度进行计算。但是,当我运行它时,有效带宽使用率突然下降到约14.4 GB / s,仅比完全双倍版本快12%。怎么能理解这一点?
我正在使用带有GCC 4.9.3的Ubuntu 14.04。
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上的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
最佳答案
必须立即进行转换的双浮点循环不能这么快地发出。通过展开一些循环,gcc可能会做得更好。
您的i7-2620M is a dual-core with hyperthreading。当瓶颈是CPU uop吞吐量,而不是分支错误预测,高速缓存未命中或什至只是较长的延迟链时,超线程无济于事。仅凭标量操作来饱和内存带宽并不容易。
从Godbolt Compiler Explorer上的代码的asm输出中可以看到:gcc 5.3具有相同的内部循环BTW,因此在这种情况下,使用旧版gcc不会造成太多损失。
double-double版本内部循环(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个融合域uops,每两个时钟可以发出一次iter。它也可以执行得很快:三个uop是负载,SnB每2个时钟可以完成4个负载。剩下5个ALU指令(因为SnB不能消除在IvB中引入的重命名阶段中的reg-reg移动)。无论如何,在单个执行端口上没有明显的瓶颈。 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个融合域uops:可以每三个时钟发出一次迭代。我认为这是唯一的瓶颈,因为它只是一个额外的需要执行端口的ALU uop。 (SnB handles zeroing-idioms in the rename stage, so
xorpd
doesn't need one)。 cvtss2sd
是2 uop指令,即使gcc使用单寄存器寻址模式,显然也不能微熔丝。每个时钟的吞吐量为一个。 (根据Agner Fog's testing的说法,在Haswell上,与寄存器src和dest一起使用时,它是一条2 uop指令,而在Skylake上,吞吐量已降低为每2个时钟1个。)这仍然不是此循环的瓶颈。不过,Skylake。在Haswell / Skylake上,它仍然是10个融合域的对象,这仍然是瓶颈。-funroll-loops
应该有助于gcc 4.9.3gcc的表现不错,代码如下
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个融合域uops,并且这不是一个微小的循环,每个第3个周期仅发出2 uops(因为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。关于c++ - 单/ double SpMV在CPU上的性能,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/36559862/