c++ - 表达式模板与手写代码

Question

我目前正在编写一个 C++ 模板表达式库，并在汇编级别将一些实例化与手写代码进行比较。手写函数如下:

spinor multiply(vector const& a, vector const& b)
{
        spinor result = {
                a.at<1>() * b.at<1>() - a.at<2>() * b.at<2>()
                          - a.at<4>() * b.at<4>() - a.at<8>() * b.at<8>(),
                a.at<1>() * b.at<2>() - a.at<2>() * b.at<1>(),
                a.at<1>() * b.at<4>() - a.at<4>() * b.at<1>(),
                a.at<1>() * b.at<8>() - a.at<8>() * b.at<1>(),
                a.at<2>() * b.at<4>() - a.at<4>() * b.at<2>(),
                a.at<2>() * b.at<8>() - a.at<8>() * b.at<2>(),
                a.at<4>() * b.at<8>() - a.at<8>() * b.at<4>()
        };

        return result;
}

vector 类只是四个 double 的包装器，可以使用 at<index>() 成员函数读取。由于设计决策，四个组件的索引是 1, 2, 4, 8 ，可以使用 at<index>() 而不是通常的 0, 1, 2, 3 访问。

此函数的目的是返回两个 vector 相乘的结果(在 Minkowski 空间中)。如果您熟悉几何代数，您将看到点积(result 的第一个组件，在 a 和 b 交换下对称)和楔积(其余组件，在 a 和 b 交换下反对称)。如果你不熟悉几何代数，就把这个函数当作 vector 乘法的处方。

如果我用 GCC 4.7 编译上面的函数并查看 objdump -SC a.out 给出的反汇编，这会给我以下输出:

400bc0: movsd  0x8(%rsi),%xmm6
400bc5: mov    %rdi,%rax
400bc8: movsd  (%rsi),%xmm8
400bcd: movsd  0x8(%rdx),%xmm5
400bd2: movapd %xmm6,%xmm9
400bd7: movsd  (%rdx),%xmm7
400bdb: movapd %xmm8,%xmm0
400be0: mulsd  %xmm5,%xmm9
400be5: movsd  0x10(%rsi),%xmm4
400bea: mulsd  %xmm7,%xmm0
400bee: movsd  0x10(%rdx),%xmm1
400bf3: movsd  0x18(%rdx),%xmm3
400bf8: movsd  0x18(%rsi),%xmm2
400bfd: subsd  %xmm9,%xmm0
400c02: movapd %xmm4,%xmm9
400c07: mulsd  %xmm1,%xmm9
400c0c: subsd  %xmm9,%xmm0
400c11: movapd %xmm3,%xmm9
400c16: mulsd  %xmm2,%xmm9
400c1b: subsd  %xmm9,%xmm0
400c20: movapd %xmm6,%xmm9
400c25: mulsd  %xmm7,%xmm9
400c2a: movsd  %xmm0,(%rdi)
400c2e: movapd %xmm5,%xmm0
400c32: mulsd  %xmm8,%xmm0
400c37: subsd  %xmm9,%xmm0
400c3c: movapd %xmm4,%xmm9
400c41: mulsd  %xmm7,%xmm9
400c46: mulsd  %xmm2,%xmm7
400c4a: movsd  %xmm0,0x8(%rdi)
400c4f: movapd %xmm1,%xmm0
400c53: mulsd  %xmm8,%xmm0
400c58: mulsd  %xmm3,%xmm8
400c5d: subsd  %xmm9,%xmm0
400c62: subsd  %xmm7,%xmm8
400c67: movapd %xmm4,%xmm7
400c6b: mulsd  %xmm5,%xmm7
400c6f: movsd  %xmm0,0x10(%rdi)
400c74: mulsd  %xmm2,%xmm5
400c78: movapd %xmm1,%xmm0
400c7c: mulsd  %xmm6,%xmm0
400c80: movsd  %xmm8,0x18(%rdi)
400c86: mulsd  %xmm3,%xmm6
400c8a: mulsd  %xmm2,%xmm1
400c8e: mulsd  %xmm4,%xmm3
400c92: subsd  %xmm7,%xmm0
400c96: subsd  %xmm5,%xmm6
400c9a: subsd  %xmm1,%xmm3
400c9e: movsd  %xmm0,0x20(%rdi)
400ca3: movsd  %xmm6,0x28(%rdi)
400ca8: movsd  %xmm3,0x30(%rdi)
400cad: retq   
400cae: nop
400caf: nop

这对我来说看起来很不错 - 第一个 ( %rsi ) 和第二个 ( %rdx ) vector 的组件仅被访问一次，实际计算仅在寄存器中完成。最后，将结果写入寄存器 %rdi 中的地址。由于这是第一个参数寄存器，我认为这里采用了返回值优化。

将此与上面函数的表达式模板版本的以下列表进行比较:

400cb0: mov    (%rsi),%rdx
400cb3: mov    0x8(%rsi),%rax
400cb7: movsd  0x1f1(%rip),%xmm4        # 400eb0 <_IO_stdin_used+0x10>
400cbe: 
400cbf: movsd  0x10(%rdx),%xmm3
400cc4: movsd  0x18(%rdx),%xmm0
400cc9: mulsd  0x10(%rax),%xmm3
400cce: xorpd  %xmm4,%xmm0
400cd2: mulsd  0x18(%rax),%xmm0
400cd7: movsd  0x8(%rdx),%xmm2
400cdc: movsd  (%rdx),%xmm1
400ce0: mulsd  0x8(%rax),%xmm2
400ce5: mulsd  (%rax),%xmm1
400ce9: subsd  %xmm3,%xmm0
400ced: subsd  %xmm2,%xmm0
400cf1: addsd  %xmm0,%xmm1
400cf5: movsd  %xmm1,(%rdi)
400cf9: movsd  (%rdx),%xmm0
400cfd: movsd  0x8(%rdx),%xmm1
400d02: mulsd  0x8(%rax),%xmm0
400d07: mulsd  (%rax),%xmm1
400d0b: subsd  %xmm1,%xmm0
400d0f: movsd  %xmm0,0x8(%rdi)
400d14: movsd  (%rdx),%xmm0
400d18: movsd  0x10(%rdx),%xmm1
400d1d: mulsd  0x10(%rax),%xmm0
400d22: mulsd  (%rax),%xmm1
400d26: subsd  %xmm1,%xmm0
400d2a: movsd  %xmm0,0x10(%rdi)
400d2f: movsd  0x8(%rdx),%xmm0
400d34: movsd  0x10(%rdx),%xmm1
400d39: mulsd  0x10(%rax),%xmm0
400d3e: mulsd  0x8(%rax),%xmm1
400d43: subsd  %xmm1,%xmm0
400d47: movsd  %xmm0,0x18(%rdi)
400d4c: movsd  (%rdx),%xmm0
400d50: movsd  0x18(%rdx),%xmm1
400d55: mulsd  0x18(%rax),%xmm0
400d5a: mulsd  (%rax),%xmm1
400d5e: subsd  %xmm1,%xmm0
400d62: movsd  %xmm0,0x20(%rdi)
400d67: movsd  0x8(%rdx),%xmm0
400d6c: movsd  0x18(%rdx),%xmm1
400d71: mulsd  0x18(%rax),%xmm0
400d76: mulsd  0x8(%rax),%xmm1
400d7b: subsd  %xmm1,%xmm0
400d7f: movsd  %xmm0,0x28(%rdi)
400d84: movsd  0x10(%rdx),%xmm0
400d89: movsd  0x18(%rdx),%xmm1
400d8e: mulsd  0x18(%rax),%xmm0
400d93: mulsd  0x10(%rax),%xmm1
400d98: subsd  %xmm1,%xmm0
400d9c: movsd  %xmm0,0x30(%rdi)
400da1: retq   

这个函数的签名是

spinor<product<vector, vector>>(product<vector, vector> const&)

我希望你相信我，两个版本都给出相同的结果。前两行提取第一个和第二个 vector ，它们作为引用存储在 product 中。我想知道以下事情:

• movsd 0x1f1(%rip),%xmm4 与 xorpd %xmm4,%xmm0 的结合有什么作用？我已经发现这叫做“RIP 相对寻址”，参见 http://www.x86-64.org/documentation/assembly.html
• 为什么 GCC 不使用更多的寄存器，例如缓存被读取四次的0x10(%rax)？

我还通过生成 100000000 个随机 vector 并花费这两个函数所需的时间来对这两个函数进行基准测试:

ET: 7.5 sec
HW: 6.8 sec

手写函数快10%左右。有没有人有使用表达式模板的经验并且知道如何使它们的性能更接近于手写模板？

Answer 1

如果我们确实知道地址 0x400eb0 的内容就很清楚了，但我怀疑它是 0x8000 0000 0000 0000 8000 0000 0000 0000，或类似的(可能带有前导 0，因为代码未矢量化)，写为 128 位 int。

在这种情况下，xorpd 会改变第二个操作数的符号。

为什么不缓存寄存器读取 - 最好在 gcc-help 邮件列表上询问这个问题。编译器可能无法证明这两个 vector 或中间结果没有混叠。

但与一般观点相反，编译器并不总是完美优化，而只是优于 90%(或 99%？)的所有程序员(如果他们尝试编写汇编)，有时(很少)他们会生成非常慢的代码.

但是您的方法非常好 - 如果您想优化，基准测试和查看生成的目标代码是正确的做法。

PS:他们的代码可能会通过使用 vector 指令(mulpd 而不是 mulsd)来加速，它会一次乘以两个或四个 double )，又名 SSE或 AVX。但是需要一些指令将值洗牌到寄存器中的正确位置，因此增益总是慢于两倍或四倍。