为什么这段代码在对循环携带的加法进行强度降低的乘法运算后执行得更慢?

我在阅读Agner Fog的优化手册时，看到了这个例子:

double data[LEN];

void compute()
{
    const double A = 1.1, B = 2.2, C = 3.3;

    int i;
    for(i=0; i<LEN; i++) {
        data[i] = A*i*i + B*i + C;
    }
}

Agner指出，有一种方法可以优化这段代码——通过实现循环可以避免使用昂贵的乘法，而是使用每次迭代应用的“增量”。

我用一张纸来证实这个理论，首先……

.．.当然，他是对的——在每次循环迭代中，我们都可以根据旧的结果计算出新的结果，通过添加一个“delta”。这个增量从值“A+B”开始，然后在每一步增加“2*A”。

所以我们将代码更新为如下所示:

void compute()
{
    const double A = 1.1, B = 2.2, C = 3.3;
    const double A2 = A+A;
    double Z = A+B;
    double Y = C;

    int i;
    for(i=0; i<LEN; i++) {
        data[i] = Y;
        Y += Z;
        Z += A2;
    }
}

就操作复杂性而言，这两个版本的函数确实存在显著差异。在我们的cpu中，与加法相比，乘法的速度要慢得多。我们已经替换了3个乘法和2个加法…只有2个补充!

所以我继续添加一个循环来执行计算很多次-然后保持执行所需的最小时间:

unsigned long long ts2ns(const struct timespec *ts)
{
    return ts->tv_sec * 1e9 + ts->tv_nsec;
}

int main(int argc, char *argv[])
{
    unsigned long long mini = 1e9;
    for (int i=0; i<1000; i++) {
        struct timespec t1, t2;
        clock_gettime(CLOCK_MONOTONIC_RAW, &t1);
        compute();
        clock_gettime(CLOCK_MONOTONIC_RAW, &t2);
        unsigned long long diff = ts2ns(&t2) - ts2ns(&t1);
        if (mini > diff) mini = diff;
    }
    printf("[-] Took: %lld ns.\n", mini);
}

我编译了两个版本，运行它们…看看这个:

gcc -O3 -o 1 ./code1.c

gcc -O3 -o 2 ./code2.c

./1

[-] Took: 405858 ns.

./2

[-] Took: 791652 ns.

这可真出乎意料。由于我们报告了最小的执行时间，因此我们丢弃了由操作系统各个部分引起的“噪音”。我们还小心地在一台完全不做任何事情的机器上运行。结果或多或少是可重复的-重新运行两个二进制文件显示这是一个一致的结果:

for i in {1..10} ; do ./1 ; done

[-] Took: 406886 ns.
[-] Took: 413798 ns.
[-] Took: 405856 ns.
[-] Took: 405848 ns.
[-] Took: 406839 ns.
[-] Took: 405841 ns.
[-] Took: 405853 ns.
[-] Took: 405844 ns.
[-] Took: 405837 ns.
[-] Took: 406854 ns.

for i in {1..10} ; do ./2 ; done

[-] Took: 791797 ns.
[-] Took: 791643 ns.
[-] Took: 791640 ns.
[-] Took: 791636 ns.
[-] Took: 791631 ns.
[-] Took: 791642 ns.
[-] Took: 791642 ns.
[-] Took: 791640 ns.
[-] Took: 791647 ns.
[-] Took: 791639 ns.

接下来要做的唯一一件事就是看看编译器为这两个版本分别创建了什么样的代码。

objdump -d - s显示了compute的第一个版本——“愚蠢的”，但以某种方式快速的代码——有一个像这样的循环:

第二个优化版本呢?它只增加了两个功能。

我不知道你们怎么想，但就我自己而言，我……困惑。第二个版本的指令大约减少了4倍，其中两个主要的指令只是基于sse的添加(addsd)。第一个版本，不仅有4倍多的指令…它还充满了(正如预期的那样)乘法(mulpd)。

我承认我没有预料到那个结果。不是因为我是阿格纳的粉丝(我是，但这无关紧要)。

你知道我错过了什么吗?我在这里犯了什么错误，可以解释速度上的差异吗?请注意，我已经在Xeon W5580和Xeon E5-1620上进行了测试-在这两个版本中，第一个(哑)版本比第二个版本快得多。

为了更容易地重现结果，这两个版本的代码有两个gist:愚蠢，但在某种程度上更快;优化，但在某种程度上更慢。

附注:请不要评论浮点精度问题;这不是问题的重点。

当前回答

如果只使用加法作为优化，就会错过(新cpu的)乘法管道的所有GFLOPS，而循环携带的依赖关系会因为停止自动向量化而使情况变得更糟。如果它是自向量化的，它将比乘法+加法快得多。并且每个数据的能源效率更高(只添加比mul + add更好)。

另一个问题是数组的末尾会收到更多的舍入错误，因为累积的加法数量。但在非常大的数组(除非数据类型变成float)之前，它不应该是可见的。

当你应用Horner的方案与GCC构建选项(在较新的cpu上)-std=c++20 -O3 -march=native -mavx2 -mprefer-vector-width=256 -ftree-vectorize -fno-math-errno，

void f(double * const __restrict__ data){
    double A=1.1,B=2.2,C=3.3;
    for(int i=0; i<1024; i++) {
        double id = double(i);
        double result =  A;
        result *=id;
        result +=B;
        result *=id;
        result += C;
        data[i] = result;
    }
}

编译器生成:

.L2:
    vmovdqa ymm0, ymm2
    vcvtdq2pd       ymm1, xmm0
    vextracti128    xmm0, ymm0, 0x1
    vmovapd ymm7, ymm1
    vcvtdq2pd       ymm0, xmm0
    vmovapd ymm6, ymm0
    vfmadd132pd     ymm7, ymm4, ymm5
    vfmadd132pd     ymm6, ymm4, ymm5
    add     rdi, 64
    vpaddd  ymm2, ymm2, ymm8
    vfmadd132pd     ymm1, ymm3, ymm7
    vfmadd132pd     ymm0, ymm3, ymm6
    vmovupd YMMWORD PTR [rdi-64], ymm1
    vmovupd YMMWORD PTR [rdi-32], ymm0
    cmp     rax, rdi
    jne     .L2
    vzeroupper
    ret

-mavx512f -mprefer-vector-width=512:

.L2:
    vmovdqa32       zmm0, zmm3
    vcvtdq2pd       zmm4, ymm0
    vextracti32x8   ymm0, zmm0, 0x1
    vcvtdq2pd       zmm0, ymm0
    vmovapd zmm2, zmm4
    vmovapd zmm1, zmm0
    vfmadd132pd     zmm2, zmm6, zmm7
    vfmadd132pd     zmm1, zmm6, zmm7
    sub     rdi, -128
    vpaddd  zmm3, zmm3, zmm8
    vfmadd132pd     zmm2, zmm5, zmm4
    vfmadd132pd     zmm0, zmm5, zmm1
    vmovupd ZMMWORD PTR [rdi-128], zmm2
    vmovupd ZMMWORD PTR [rdi-64], zmm0
    cmp     rax, rdi
    jne     .L2
    vzeroupper
    ret

由于mul+add连接到单个FMA，所有浮点运算都是“打包”的向量形式和更少的指令(它是一个两次展开的版本)。每64字节数据16条指令(AVX-512则为128字节)。

Horner方案的另一个优点是，它在FMA指令中计算精度更高，每次循环迭代只有O(1)个操作，因此对于较长的数组，它不会累积那么多错误。

我认为Agner Fog的优化手册中的优化一定来自《雷神之锤III》的快速平方根近似，这在x87中是有意义的，但在SSE中被淘汰了。在那个时候，SIMD不够宽，不能产生太大的差异，而且对于sqrt()和除法来说也很慢，但是有rsqrtss。手册上写着2004年的版权，所以赛扬的cpu是SSE而不是FMA。第一个AVX桌面CPU的推出要晚得多，FMA甚至比它还要晚。

下面是另一个强度降低的版本(用于id值):

void f(double * const __restrict__ data){

    double B[]={2.2,2.2,2.2,2.2,2.2,2.2,2.2,2.2,
    2.2,2.2,2.2,2.2,2.2,2.2,2.2,2.2};
    double C[]={3.3,3.3,3.3,3.3,3.3,3.3,3.3,3.3,
    3.3,3.3,3.3,3.3,3.3,3.3,3.3,3.3};
    double id[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};
    for(long long i=0; i<1024; i+=16) {
        double result[]={1.1,1.1,1.1,1.1,1.1,1.1,1.1,1.1,
                        1.1,1.1,1.1,1.1,1.1,1.1,1.1,1.1};
        // The same thing, just with explicit auto-vectorization help
        for(int j=0;j<16;j++)
        {
            result[j] *=id[j];
            result[j] +=B[j];
            result[j] *=id[j];
            result[j] += C[j];
            data[i+j] = result[j];
        }

        // strength reduction
        for(int j=0;j<16;j++)
        {
            id[j] += 16.0;
        }
    }
}

组装:

.L2:
    vmovapd zmm3, zmm0
    vmovapd zmm2, zmm1
    sub     rax, -128
    vfmadd132pd     zmm3, zmm6, zmm7
    vfmadd132pd     zmm2, zmm6, zmm7
    vfmadd132pd     zmm3, zmm5, zmm0
    vfmadd132pd     zmm2, zmm5, zmm1
    vaddpd  zmm0, zmm0, zmm4
    vaddpd  zmm1, zmm1, zmm4
    vmovupd ZMMWORD PTR [rax-128], zmm3
    vmovupd ZMMWORD PTR [rax-64], zmm2
    cmp     rdx, rax
    jne     .L2
    vzeroupper
    ret

当数据、A、B和C数组通过alignas(64)和足够小的数据数组大小进行对齐时，它以每个元素速度0.26个周期运行。

2022-05-20 13:06:39

其他回答

如果你需要这段代码来快速运行，或者如果你很好奇，你可以尝试以下方法:

您将a[i] = f(i)的计算更改为两个加法。将其修改为使用两个加法计算a[4i] = f(4i)，使用两个加法计算a[4i+1] = f(4i+1)，依此类推。现在你有四个计算可以并行完成。

编译器很有可能会执行相同的循环展开和向量化，并且您有相同的延迟，但对于四个操作，而不是一个。

2022-05-20 11:48:24

The main difference here is loop dependencies. The loop in the second case is dependent—operations in the loop depend on the previous iteration. This means that each iteration can't even start until the previous iteration finishes. In the first case, the loop body is fully independent—everything in the loop body is self contained, depending solely on the iteration counter and constant values. This means that the loop can be computed in parallel—multiple iterations can all work at the same time. This then allows the loop to be trivially unrolled and vectorized, overlapping many instructions.

如果您要查看性能计数器(例如，perf stat ./1)，您将看到第一个循环除了运行得更快之外，每个周期(IPC)还运行了更多的指令。相比之下，第二个循环有更多的依赖周期——当CPU无所事事，等待指令完成，然后才能发出更多指令时。

第一个可能会导致内存带宽瓶颈，特别是如果你让编译器在Sandy Bridge (gcc -O3 -march=native)上使用AVX自动向量化，如果它设法使用256位向量的话。此时，IPC将会中断，特别是对于一个对于L3缓存来说太大的输出数组。

注意:展开和向量化不需要独立的循环——当存在(某些)循环依赖关系时可以执行它们。然而，这样做难度更大，回报也更少。因此，如果你想看到从向量化得到的最大加速，它有助于在可能的情况下删除循环依赖。

2022-05-19 16:42:51

理解您所看到的性能差异的关键在于向量化。是的，基于加法的解决方案在其内部循环中只有两条指令，但重要的区别不在于循环中有多少条指令，而在于每条指令执行了多少工作。

在第一个版本中，输出完全依赖于输入:每个数据[i]只是i本身的函数，这意味着每个数据[i]可以以任何顺序计算:编译器可以向前、向后、横向，无论如何，你仍然会得到相同的结果——除非你从另一个线程观察内存，否则你永远不会注意到数据是如何被处理的。

在第二个版本中，输出不依赖于i——它依赖于最后一次循环中的A和Z。

如果我们将这些循环的主体表示为小的数学函数，它们将有非常不同的整体形式:

f(i) - f(Y, Z) -> (di, Y, Z)

在后一种形式中，没有实际依赖于i -你可以计算函数值的唯一方法是通过从函数的最后一次调用中知道前一个Y和Z，这意味着函数形成了一个链-你不能做下一个，直到你做了前一个。

Why does that matter? Because the CPU has vector parallel instructions that each can perform two, four, or even eight arithmetic operations at the same time! (AVX CPUs can do even more in parallel.) That's four multiplies, four adds, four subtracts, four comparisons — four whatevers! So if the output you're trying to compute is only dependent on the input, then you can safely do two, four, or even eight at a time — it doesn't matter if they're forward or backward, since the result is the same. But if the output is dependent on previous computation, then you're stuck doing it in serial form — one at a time.

这就是为什么“较长的”代码在性能方面更胜一筹。尽管它有更多的设置，而且实际上做了更多的工作，但大部分工作是并行完成的:它不仅仅是在循环的每次迭代中计算数据[i]，而是同时计算数据[i]，数据[i+1]，数据[i+2]和数据[i+3]，然后跳转到下一个4个数据的集合。

为了扩展一下我在这里的意思，编译器首先将原始代码转换成这样:

int i;
for (i = 0; i < LEN; i += 4) {
    data[i+0] = A*(i+0)*(i+0) + B*(i+0) + C;
    data[i+1] = A*(i+1)*(i+1) + B*(i+1) + C;
    data[i+2] = A*(i+2)*(i+2) + B*(i+2) + C;
    data[i+3] = A*(i+3)*(i+3) + B*(i+3) + C;
}

你可以说服自己，如果你眯着眼睛看，它会和原来的一样。它之所以能做到这一点，是因为所有这些相同的垂直操作符行:所有这些*和+操作都是相同的操作，只是在不同的数据上执行——而且CPU有特殊的内置指令，可以在同一时间对不同的数据执行多个*或多个+操作，每个操作只需一个时钟周期。

注意快速解决方案addpd和mulpd指令中的字母p，以及慢速解决方案addsd指令中的字母s。这是“添加打包双打”和“增加打包双打”，而不是“添加单个双打”。

不仅如此，它看起来像编译器也部分展开了循环-循环不是每次迭代只执行两个值，而是实际上执行了四个值，并且交错操作以避免依赖关系和停顿，所有这些都减少了汇编代码测试i < 1000的次数。

All of this only works, though, if there are no dependencies between iterations of the loop: If the only thing that determines what happens for each data[i] is i itself. If there are dependencies, if data from the last iteration influences the next one, then the compiler may be so constrained by them that it can't alter the code at all — instead of the compiler being able to use fancy parallel instructions or clever optimizations (CSE, strength reduction, loop unrolling, reordering, et al.), you get out code that's exactly what you put in — add Y, then add Z, then repeat.

但是这里，在代码的第一个版本中，编译器正确地识别出数据中没有依赖关系，并发现它可以并行地完成工作，所以它做到了，这就是所有不同的地方。

2022-05-19 16:29:59

通过手动将代码并行化，你似乎可以鱼与熊掌兼得:

double A4 = A+A+A+A;
double Z = 3A+B;
double Y1 = C;
double Y2 = A+B+C;

int i;
// ... setup unroll when LEN is odd...

for(i=0; i<LEN; i++) {
    data[i] = Y1;
    data[++i] = Y2;
    Y1 += Z;
    Y2 += Z;
    Z += A4;
}

它可能不完全像所写的那样具有功能性，但您可以理解它的思想:展开循环，以便每个依赖于数据的路径都可以并行执行。对于所考虑的机器，四步展开应该能够实现最大的性能，当然，在软件中对体系结构进行硬编码时，您可以获得所有有趣的东西。

2022-05-29 19:02:55