获取 std::complex 的组件作为参考

Question

std::real(my_complex) 和 my_complex.real() 是否复制真实部分？有没有一种方法可以通过引用而不是值来访问？

对于背景，我正在编写一些性能关键代码。在紧密的循环中，我必须做一些复杂的 * 实数乘法。我发现执行两个实数乘法比执行复数乘法更快，因为我知道其中一个操作数是实数。为了支持实数乘法，我将复数数据存储为 SOA、std::complex<std::vector<short>>。也许这是个坏主意，但我认为这会让 reader 明白这是存储为数组结构的复杂数据。

无论如何，在紧密循环中，我会执行以下操作：

std::real(complex_data)[0] * all_real_data[0]
std::imag(complex_data)[0] * all_real_data[0]

事实证明，在 CPU 使用报告中，真实和图像查找是大问题。

我试了complex_data.real()[0] * all_real_data[0]，好像也没什么区别

然后我从循环中提取了 real/imag 尊重，比如 std::vector<short>& my_complex_real = std::real(complex_data) 速度快 2 倍。

我想子问题是 "Is SOA inside a std::complex a bad idea?"

Answer 1

std::real 和 std::complex::real 都按值给你真实的部分，这意味着他们制作了一个副本。

访问 std::complex<T> 的实部和虚部的唯一方法是将其转换为数组。如果你有

std::complex<T> foo;

然后

reinterpret_cast<T(&)[2]>(foo)[0]

给你一个真实部分的参考并且

reinterpret_cast<T(&)[2]>(foo)[1]

给你虚部参考。这是按照标准 ([complex.numbers]/4) 强制执行的，因此它不是未定义的行为。

您还应注意 std::complex 仅针对 std::complex<float>、std::complex<double> 和 std::complex<long double> 定义。 [complex.numbers]/2

未指定任何其他实例

Answer 2

我认为这里的 SOA 想法不会特别有效。我假设您必须为 std::vector 添加全局算术重载才能完成这项工作。但在内部，这也意味着有两个可调整大小的向量和两个额外的指针，对于 SOA 与 AOS 很重要的应用程序来说，这是相当多的开销。这也给出了提取实部成本高昂的原因：向量本身几乎肯定会被复制。

@NathanOliver 上面的回答提供了一种方法来获取指向 std::complex 的指针作为数组，这可能会保存复制，但我希望您至少要使用自定义 class 而不是 std::vector<short>。现实中复杂的算术很容易实现，自己做那部分可能会更快。

（Daniel H 的回答比我的要好，因为它表明规范不允许它并专门调用缓存位置。你真的不想这样做。）

Answer 3

使用 std::complex<std::vector<short>> 是未指定的行为。除非您有编译器扩展，否则 only allowed specializations 是特化 std::complex<float>、std::complex<double> 和 std::complex<long double>。其他算术类型，如 std::complex<short>，即使在理论上没有任何更强的要求，但在实践中至少更有可能获得合理的结果。

由于缓存局部性，我希望 std::vector<std::complex<short>> 会有更好的性能，即使这两种类型恰好在您的实现中都能很好地工作。

无论哪种方式，正如 NathanOliver 上面指出的那样，reinterpret_cast<T(&)[2]>(z)[0] 和 reinterpret_cast<T(&)[2]>(z)[1] 应该给出对实部和虚部的引用，但请注意，复数定义了一个 operator* 用于乘以真实类型，所以这不是必需的。

Answer 4

所以我在 godbolt 上写了这个小例子。

我很难看到复制的问题：

#include <complex>
#include <array>

std::complex<double> foo(std::complex<double> (& complex_data)[10], double (&all_real_data)[10], int i)
{
    return std::complex<double>(std::real(complex_data[i]) * all_real_data[i], 
        std::imag(complex_data[i]) * all_real_data[i]);
}

std::array<std::complex<double>, 10>
calc(std::complex<double> (& complex_data)[10], double (&all_real_data)[10])
{
    std::array<std::complex<double>, 10> result;
    for (int i = 0 ; i < 10 ; ++i)
    {
        result[i] = foo(complex_data, all_real_data, i);
    }
    return result;
}

在 gcc 上使用 -O3 编译计算的产量：

calc(std::complex<double> (&) [10], double (&) [10]):
  movsd xmm0, QWORD PTR [rdx]
  mov rax, rdi
  movsd xmm1, QWORD PTR [rsi+8]
  mulsd xmm1, xmm0
  mulsd xmm0, QWORD PTR [rsi]
  movsd QWORD PTR [rdi+8], xmm1
  movsd xmm1, QWORD PTR [rsi+24]
  movsd QWORD PTR [rdi], xmm0
  movsd xmm0, QWORD PTR [rdx+8]
  mulsd xmm1, xmm0
  mulsd xmm0, QWORD PTR [rsi+16]
  movsd QWORD PTR [rdi+24], xmm1
  movsd xmm1, QWORD PTR [rsi+40]
  movsd QWORD PTR [rdi+16], xmm0
  movsd xmm0, QWORD PTR [rdx+16]
  mulsd xmm1, xmm0
  mulsd xmm0, QWORD PTR [rsi+32]
  movsd QWORD PTR [rdi+40], xmm1
  movsd xmm1, QWORD PTR [rsi+56]
  movsd QWORD PTR [rdi+32], xmm0
  movsd xmm0, QWORD PTR [rdx+24]
  mulsd xmm1, xmm0
  mulsd xmm0, QWORD PTR [rsi+48]
  movsd QWORD PTR [rdi+56], xmm1
  movsd xmm1, QWORD PTR [rsi+72]
  movsd QWORD PTR [rdi+48], xmm0
  movsd xmm0, QWORD PTR [rdx+32]
  mulsd xmm1, xmm0
  mulsd xmm0, QWORD PTR [rsi+64]
  movsd QWORD PTR [rdi+72], xmm1
  movsd xmm1, QWORD PTR [rsi+88]
  movsd QWORD PTR [rdi+64], xmm0
  movsd xmm0, QWORD PTR [rdx+40]
  mulsd xmm1, xmm0
  mulsd xmm0, QWORD PTR [rsi+80]
  movsd QWORD PTR [rdi+88], xmm1
  movsd xmm1, QWORD PTR [rsi+104]
  movsd QWORD PTR [rdi+80], xmm0
  movsd xmm0, QWORD PTR [rdx+48]
  mulsd xmm1, xmm0
  mulsd xmm0, QWORD PTR [rsi+96]
  movsd QWORD PTR [rdi+104], xmm1
  movsd xmm1, QWORD PTR [rsi+120]
  movsd QWORD PTR [rdi+96], xmm0
  movsd xmm0, QWORD PTR [rdx+56]
  mulsd xmm1, xmm0
  mulsd xmm0, QWORD PTR [rsi+112]
  movsd QWORD PTR [rdi+120], xmm1
  movsd xmm1, QWORD PTR [rsi+136]
  movsd QWORD PTR [rdi+112], xmm0
  movsd xmm0, QWORD PTR [rdx+64]
  mulsd xmm1, xmm0
  mulsd xmm0, QWORD PTR [rsi+128]
  movsd QWORD PTR [rdi+136], xmm1
  movsd xmm1, QWORD PTR [rsi+152]
  movsd QWORD PTR [rdi+128], xmm0
  movsd xmm0, QWORD PTR [rdx+72]
  mulsd xmm1, xmm0
  mulsd xmm0, QWORD PTR [rsi+144]
  movsd QWORD PTR [rdi+152], xmm1
  movsd QWORD PTR [rdi+144], xmm0
  ret

使用 -march=native 我们接触内存的次数更少

calc(std::complex<double> (&) [10], double (&) [10]):
  vmovupd ymm1, YMMWORD PTR [rsi]
  vmovupd ymm0, YMMWORD PTR [rsi+32]
  mov rax, rdi
  vmovupd ymm3, YMMWORD PTR [rdx]
  vunpckhpd ymm2, ymm1, ymm0
  vunpcklpd ymm0, ymm1, ymm0
  vpermpd ymm2, ymm2, 216
  vpermpd ymm0, ymm0, 216
  vmulpd ymm0, ymm0, ymm3
  vmulpd ymm2, ymm2, ymm3
  vpermpd ymm1, ymm0, 68
  vpermpd ymm0, ymm0, 238
  vpermpd ymm3, ymm2, 68
  vpermpd ymm2, ymm2, 238
  vshufpd ymm1, ymm1, ymm3, 12
  vshufpd ymm0, ymm0, ymm2, 12
  vmovupd YMMWORD PTR [rdi], ymm1
  vmovupd ymm1, YMMWORD PTR [rsi+64]
  vmovupd YMMWORD PTR [rdi+32], ymm0
  vmovupd ymm0, YMMWORD PTR [rsi+96]
  vmovupd ymm3, YMMWORD PTR [rdx+32]
  vunpckhpd ymm2, ymm1, ymm0
  vunpcklpd ymm0, ymm1, ymm0
  vpermpd ymm2, ymm2, 216
  vpermpd ymm0, ymm0, 216
  vmulpd ymm0, ymm0, ymm3
  vmulpd ymm2, ymm2, ymm3
  vpermpd ymm1, ymm0, 68
  vpermpd ymm0, ymm0, 238
  vpermpd ymm3, ymm2, 68
  vpermpd ymm2, ymm2, 238
  vshufpd ymm1, ymm1, ymm3, 12
  vshufpd ymm0, ymm0, ymm2, 12
  vmovupd YMMWORD PTR [rdi+64], ymm1
  vmovupd YMMWORD PTR [rdi+96], ymm0
  vmovsd xmm0, QWORD PTR [rdx+64]
  vmulsd xmm1, xmm0, QWORD PTR [rsi+136]
  vmulsd xmm0, xmm0, QWORD PTR [rsi+128]
  vmovsd QWORD PTR [rdi+136], xmm1
  vmovsd QWORD PTR [rdi+128], xmm0
  vmovsd xmm0, QWORD PTR [rdx+72]
  vmulsd xmm1, xmm0, QWORD PTR [rsi+152]
  vmulsd xmm0, xmm0, QWORD PTR [rsi+144]
  vmovsd QWORD PTR [rdi+152], xmm1
  vmovsd QWORD PTR [rdi+144], xmm0
  vzeroupper
  ret

获取 std::complex 的组件作为参考

Fetch component of std::complex as reference

c++

std