g++与手动优化在复数乘法中的比较

9
在我们的代码库中,有许多类似于j*ω*X的操作,其中j是虚数单位,ω是实数,X是复数。实际上,很多循环都可能看起来像这样:
#include <complex>
#include <vector>

void mult_jomega(std::vector<std::complex<double> > &vec, double omega){
    std::complex<double> jomega(0.0, omega);
    for (auto &x : vec){
        x*=jomega;
    }
}

然而,我们利用实部为零的事实,将乘法写成以下形式:

jomega

void mult_jomega_smart(cvector &vec, double omega){
    for (auto &x : vec){
        x={-omega*x.imag(), omega*x.real()};
    }
}

一开始我对这个“智能”版本持怀疑态度,因为:

  1. 它更难理解。
  2. 出错的可能性更大。
  3. “编译器会进行优化。”

然而,由于一些性能回归测试表明,第三个论据并不成立。在比较这两个函数时(请参见下面的代码),使用-O2-O3编译选项时,“智能”版本始终表现出色:

size    orig(musec)   smart(musec)  speedup
10      0.039928      0.0117551     3.39665
100     0.328564      0.0861379     3.81439
500     1.62269       0.417475      3.8869
1000    3.33012       0.760515      4.37877
2000    6.46696       1.56048       4.14422
10000   32.2827       9.2361        3.49528
100000  326.828       115.158       2.8381
500000  1660.43       850.415       1.95249

智能版本在我的机器上(gcc-5.4)大约快了4倍,只有当任务随着数组大小的增加变得越来越受内存限制时,加速才会下降到2倍左右。
我的问题是,为什么编译器不能优化不太智能但更易读的版本,毕竟编译器可以看到jomega的实部为零?是否可以通过提供一些额外的编译标志来帮助编译器进行优化?
注:其他编译器也存在加速。
compiler      speedup
g++-5.4          4
g++-7.2          4
clang++-3.8      2  [original version 2-times faster than gcc]

代码清单:

mult.cpp - 防止内联:

#include <complex>
#include <vector>

typedef std::vector<std::complex<double> > cvector;
void mult_jomega(cvector &vec, double omega){
    std::complex<double> jomega(0.0, omega);
    for (auto &x : vec){
        x*=jomega;
    }
}

void mult_jomega_smart(cvector &vec, double omega){
    for (auto &x : vec){
        x={-omega*x.imag(), omega*x.real()};
    }
}

main.cpp:

#include <chrono>
#include <complex>
#include <vector>
#include <iostream>

typedef std::vector<std::complex<double> > cvector;
void mult_jomega(cvector &vec, double omega);
void mult_jomega2(cvector &vec, double omega);
void mult_jomega_smart(cvector &vec, double omega);


const size_t N=100000;   //10**5
const double OMEGA=1.0;//use 1, so nothing changes -> no problems with inf & Co

void compare_results(const cvector &vec){
   cvector m=vec;
   cvector m_smart=vec;
   mult_jomega(m, 5.0);
   mult_jomega_smart(m_smart,5.0);
   std::cout<<m[0]<<" vs "<<m_smart[0]<<"\n";
   std::cout<< (m==m_smart ? "equal!" : "not equal!")<<"\n";

}

void test(size_t vector_size){

     cvector vec(vector_size, std::complex<double>{1.0, 1.0});

     //compare results, triger if in doubt
     //compare_results(vec);


     //warm_up, just in case:
     for(size_t i=0;i<N;i++)
        mult_jomega(vec, OMEGA);

     //test mult_jomega:
     auto begin = std::chrono::high_resolution_clock::now();
     for(size_t i=0;i<N;i++)
        mult_jomega(vec, OMEGA);
     auto end = std::chrono::high_resolution_clock::now();
     auto time_jomega=std::chrono::duration_cast<std::chrono::nanoseconds>(end-begin).count()/1e3;


     //test mult_jomega_smart:
     begin = std::chrono::high_resolution_clock::now();
     for(size_t i=0;i<N;i++)
        mult_jomega_smart(vec, OMEGA);
     end = std::chrono::high_resolution_clock::now();
     auto time_jomega_smart=std::chrono::duration_cast<std::chrono::nanoseconds>(end-begin).count()/1e3;

     double speedup=time_jomega/time_jomega_smart;
     std::cout<<vector_size<<"\t"<<time_jomega/N<<"\t"<<time_jomega_smart/N<<"\t"<<speedup<<"\n";
}


int main(){
   std::cout<<"N\tmult_jomega(musec)\tmult_jomega_smart(musec)\tspeedup\n";    
   for(const auto &size : std::vector<size_t>{10,100,500,1000,2000,10000,100000,500000})
        test(size);          
}

构建和运行:

g++ main.cpp mult.cpp -O3 -std=c++11 -o mult_test
./mult_test

1
jomega 不是 const。此外,operator * 很可能采用 const & 并允许 const_cast 去除 const 并更改 jomega,因此不能确定 jomega 的实部是否为 0.0。尝试将 const 添加到 jomega 中,看看是否有任何变化。 - nwp
2
添加 const 没有太大作用。 - Aziz
1
@nwp 我认为将某些东西声明为const并不能帮助编译器优化你的代码。0.0是一个字面量,常量传播应该无论如何都会发生。有一个旧的gotw声明说const对于优化没有任何作用(当然constexpr有):http://www.gotw.ca/gotw/081.htm。 - Jens
@Jens 通常来说,const 对于优化并没有什么作用(比如,一个函数使用 const & 或者 & 并没有太大区别),但是顶层的 const 实际上是意味着不可变的,因为强制转换会导致未定义行为。话虽如此,显然我错了,在这里也没有什么作用。 - nwp
3个回答

8

使用标志-ffast-math进行编译可以获得更快的性能。

N       mult_jomega(musec)      mult_jomega_smart(musec)        speedup
10      0.00860809              0.00818644                      1.05151
100     0.0706683               0.0693907                       1.01841
500     0.29569                 0.297323                        0.994509
1000    0.582059                0.57622                         1.01013
2000    1.30809                 1.24758                         1.0485
10000   7.37559                 7.4854                          0.98533
编辑:更具体地说,它是-funsafe-math-optimizations编译器标志。根据文档,此标志用于允许针对浮点算术的优化,假设参数和结果是有效的,并且可能违反IEEE或ANSI标准。 编辑2:更加具体地说,它是-fno-signed-zeros选项,其功能为:允许针对浮点算术进行优化,忽略零的符号。 IEEE算术规定了不同的+0.0-0.0值的行为,这就禁止简化诸如x+0.00.0*x之类的表达式(即使使用-ffinite-math-only)。此选项意味着零结果的符号不重要。

1
“-fno-signed-zeros”是一个非常好的发现。既然它已经足够了,我鼓励只使用该标志而不是“-ffast-math”,除非广泛测试表明它不会破坏任何东西。 - nwp
我理解得对吗,其中一个问题是情况 omega=0.0 可能会导致结果中的零有不同的符号,例如答案的实部符号取决于输入的实部和虚部? - ead
@ead: 正确。如果没有该标志,编译器无法优化乘以零的操作,因为 x*0.0 的结果可能是两个“不同”的数字之一:0.0-0.0。该标志打破了 IEEE 754 标准规则,并告诉编译器将零视为无符号数。因此,x*0.0 总是 0.0,编译器可以从那里进行优化。 - Aziz
思考这个问题,似乎根本问题在于涉及纯实数和复数的算术运算以及涉及纯虚数和复数的算术运算之间存在不对称性。对于浮点/复合操作,浮点数的虚部在概念上是一个没有符号的零,对于虚数/复合操作,虚数的实部的零被认为具有符号。因此,我确信您指出了常见复杂实现中的漏洞:它们应该将纯虚数实现为独立类型。 - Oliv
@ead:clang没有-fno-signed-zeros标志。 - Aziz
显示剩余3条评论

1
我对 Aziz 的回答中提到的编译器选项进行了更多调查,使用 godbolt 编译器探索器。示例代码实现了内部循环的三个版本:
  1. mult_jomega 示例。
  2. 手写版本的相同循环,其中调用 operator*= 已被替换为计算。
  3. mult_jomega_smart 示例。

来自 godbolt 的代码

// 1. mult_jomega
std::complex<double> const jomega(0.0, omega);
for (auto &x : v){
    x*=jomega;
}

// 2. hand-written mult_jomega
for (auto &x : v3){
    double x1 = x.real() * jomega.real();
    double x2 = x.imag() * jomega.imag();
    double x3 = x.real() * jomega.imag();
    double x4 = x.imag() * jomega.real();
    x = { x1 - x2 , x3 + x4};
}

// 3. mult_jomega_smart
for (auto &x : v2){
    x={-omega*x.imag(), omega*x.real()};
}

检查三个循环的汇编代码:

mult_jomega

 cmp %r13,%r12
 je 4008ac <main+0x10c>
 mov %r12,%rbx
 nopl 0x0(%rax)
 pxor %xmm0,%xmm0
 add $0x10,%rbx
 movsd -0x8(%rbx),%xmm3
 movsd -0x10(%rbx),%xmm2
 movsd 0x8(%rsp),%xmm1
 callq 400740 <__muldc3@plt>
 movsd %xmm0,-0x10(%rbx)
 movsd %xmm1,-0x8(%rbx)
 cmp %rbx,%r13
 jne 400880 <main+0xe0>

手写乘法
 cmp %rdx,%rdi
 je 40090c <main+0x16c>
 pxor %xmm3,%xmm3
 mov %rdi,%rax
 movsd 0x8(%rsp),%xmm5
 nopl 0x0(%rax,%rax,1)
 movsd (%rax),%xmm0
 movapd %xmm5,%xmm4
 movsd 0x8(%rax),%xmm1
 add $0x10,%rax
 movapd %xmm0,%xmm2
 mulsd %xmm5,%xmm0
 mulsd %xmm1,%xmm4
 mulsd %xmm3,%xmm2
 mulsd %xmm3,%xmm1
 subsd %xmm4,%xmm2
 addsd %xmm1,%xmm0
 movsd %xmm2,-0x10(%rax)
 movsd %xmm0,-0x8(%rax)
 cmp %rax,%rdx
 jne 4008d0 <main+0x130>

mult_jomega_smart

cmp %rcx,%rdx
 je 400957 <main+0x1b7>
 movsd 0x8(%rsp),%xmm2
 mov %rcx,%rax
 xorpd 0x514(%rip),%xmm2 # 400e40 <_IO_stdin_used+0x10>
 nopl 0x0(%rax)
 add $0x10,%rax
 movsd 0x8(%rsp),%xmm0
 movsd -0x8(%rax),%xmm1
 mulsd -0x10(%rax),%xmm0
 mulsd %xmm2,%xmm1
 movsd %xmm1,-0x10(%rax)
 movsd %xmm0,-0x8(%rax)
 cmp %rax,%rdx
 jne 400930 <main+0x190>

我对汇编代码的理解很有限,但是我发现:

  • operator*=没有在mult_jomega中内联
  • x1x4被计算了,尽管它们始终为0.0,因为jomega.real()==0.0

我不知道为什么operator*=没有被内联。源代码很简单,只包含三行。

当你考虑到0.0 * x == 0.0不总是对于double类型的值成立时,可以解释计算x1x4的原因。除了其他答案中提到的有符号零定义外,还有无穷大的值naninf,其中x * 0.0 = 0.0不成立。

如果使用-fno-signed-zeros-ffinite-math-only编译,优化将应用并移除对x1x4的计算。

0

正如其他答案所指出的那样,我的错误在于假设纯虚数 j*ω 与复数 0.0+j*ω 具有相同的行为方式 - 就像纯实数 1.0 与复数 1.0+0.0j 没有相同的行为方式一样,例如(使用gcc编译)。

1.0*(inf+0.0j) = inf+0.0j
(1.0 +0.0j)*(inf+0.0j) = inf - nanj

这是因为复数乘法比学校公式建议的更加复杂

基本上,C++编译器在处理复数时存在一种不对称性:虽然有纯实数(即double),但没有纯虚数。

C++标准没有定义复数乘法必须发生的方式。大多数编译器都借鉴了它们实现C99的方式。C99是第一个定义了复数操作如何进行的C格式,在Annex G中定义了它。但是,正如G.1.1中所述,它是可选的:

G.1.1...尽管这些规范经过精心设计, 但很少有现有的实践可以验证设计决策。 因此,这些规范不是规范性的,而应该被视为推荐做法...

C99还定义了纯虚数数据类型 float _Imaginary double _Imaginary long double _Imaginary (G.2),这正是我们需要的j*ω。 G.5.1定义了纯虚数的乘法语义。

xj*(v+wj) = -xw+(xv)j
(v+wj)*xj = -wx+(vx)j

即学校公式已经足够好了(与两个复数的乘法不同)。

问题在于,到目前为止,所知的编译器(gcc-9.2、clang-9.0)都没有支持_Imaginary类型的(因为它是可选的)。

因此,我的解决方案是实现一个pure_imaginary类型并按照G.5.1重载运算符。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接