为什么在强度降低乘法和循环进位加法之后，这段代码的执行速度会变慢？
tts*_*ras 320 optimization assembly x86-64 simd cpu-architecture
我正在阅读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 指出，有一种方法可以优化此代码 - 通过认识到循环可以避免使用昂贵的乘法，而是使用每次迭代应用的“增量”。
我用一张纸来证实这个理论，首先......
...当然，他是对的 - 在每次循环迭代中，我们可以通过添加“增量”，基于旧结果计算新结果。该增量从值“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 中以显着慢着称。我们用 2 个加法替换了 3 个乘法和 2 个加法......
所以我继续添加一个循环来执行compute很多次 - 然后保留执行所需的最短时间：
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上进行了测试- 在这两个版本中，第一个（哑）版本比第二个版本快得多。
为了轻松重现结果，两个版本的代码有两个要点：愚蠢但速度更快，优化但速度较慢。
PS 请不要评论浮点精度问题；这不是这个问题的重点。
Sea*_*ema 314
理解您所看到的性能差异的关键在于矢量化。是的，基于加法的解决方案在其内部循环中只有两条指令，但重要的区别不在于循环中有多少指令，而在于每条指令执行多少工作。
\n在第一个版本中，输出完全依赖于输入：每个data[i]都是其自身的函数i，这意味着每个都data[i]可以按任何顺序计算：编译器可以向前、向后、横向等等进行计算，而你\'仍然会得到相同的结果 \xe2\x80\x94 除非您从另一个线程观察该内存，否则您永远不会注意到数据的处理方式。
在第二个版本中，输出不依赖于i\xe2\x80\x94 它依赖于上次循环的A和。Z
如果我们将这些循环体表示为小数学函数，它们将具有非常不同的整体形式：
\n- \n
- f(i) -> di \n
- f(Y, Z) -> (di, Y\', Z\') \n
在后一种形式中，对 \xe2\x80\x94 没有实际依赖，i计算函数值的唯一方法是了解函数的前一个Y和Z最后一个调用，这意味着函数形成一条链\xe2\x80\x94，在完成前一项之前，您无法进行下一项。
为什么这很重要？因为CPU有向量并行指令，每个指令可以同时执行两个、四个甚至八个算术运算！（AVX CPU 可以并行执行更多操作。）这是四个乘法、四个加法、四个减法、四个比较 \xe2\x80\x94 四个无论如何！因此，如果您尝试计算的输出仅依赖于输入，那么您可以安全地一次执行两个、四个甚至八个 \xe2\x80\x94，无论它们是否是向前或向后，因为结果是相同的。但是，如果输出依赖于先前的计算，那么您将不得不以串行形式 \xe2\x80\x94 一次执行一个。
\n这就是为什么“较长”的代码在性能上获胜的原因。尽管它有更多的设置，并且实际上做了更多的工作，但大部分工作都是并行完成的：它不仅仅data[i]在循环的每次迭代中进行计算 \xe2\x80\x94它同时计算data[i]、data[i+1]、data[i+2]、 和，然后跳到下一组四个。data[i+3]
为了扩展一下我在这里的意思，编译器首先将原始代码转换为如下所示：
\nint i;\nfor (i = 0; i < LEN; i += 4) {\n data[i+0] = A*(i+0)*(i+0) + B*(i+0) + C;\n data[i+1] = A*(i+1)*(i+1) + B*(i+1) + C;\n data[i+2] = A*(i+2)*(i+2) + B*(i+2) + C;\n data[i+3] = A*(i+3)*(i+3) + B*(i+3) + C;\n}\n如果你眯着眼睛看的话，你可以说服自己，它会做与原版相同的事情。之所以这样做，是因为所有这些相同的垂直操作符：所有这些*和+操作都是相同的操作，只是在不同的数据 \xe2\x80\x94 上执行，并且 CPU 具有特殊的内置指令，可以执行多个*或+同时对不同数据进行多个操作，每个操作仅在一个时钟周期内进行。
p请注意较快解决方案 \xe2\x80\x94addpd和mulpd\xe2\x80\x94中s的说明中的字母以及较慢解决方案 \xe2\x80\x94 中的说明中的字母addsd。这就是“添加压缩双精度数”和“乘以压缩双精度数”与“添加单双精度数”。
不仅如此，编译器看起来也部分展开了循环 \xe2\x80\x94 循环每次迭代不仅执行两个值，而且实际上执行四个值，并交错操作以避免依赖性和停顿，所有这些i < 1000还减少了汇编代码必须测试的次数。
然而，所有这些只有在循环迭代之间不存在依赖关系的情况data[i]下才有效：如果唯一决定每次发生的事情是i它自己。如果存在依赖性，如果上一次迭代的数据影响下一次迭代，那么编译器可能会受到它们的限制，以至于根本无法更改代码 \xe2\x80\x94 而不是编译器能够使用花哨的并行指令或巧妙的优化（CSE、强度降低、循环展开、重新排序等），您得到的代码正是您在 \xe2\x80\x94 中放入的代码添加 Y，然后添加 Z，然后重复。
但是在这里，在代码的第一个版本中，编译器正确地识别出数据中没有依赖性，并发现它可以并行地完成工作，而且它确实做到了，这就是使所有不同之处。
\n- 这不仅仅是矢量化，还有数据依赖性。由于迭代中的延迟瓶颈，“优化”版本中的标量代码无法全速运行。这与阻止其矢量化的原因相同，但我会通过说关键是循环携带的依赖项来开始回答。缺乏这样的允许跨迭代的向量化和指令级并行性。（整数“i++”是一个循环携带的 dep，但编译器可以使用它，因为整数数学是关联的，这与没有“-ffast-math”的 FP 不同） (44认同)
- @PeterCordes我真的想在这个答案中关注“并行与串行计算”的高级概念，因为这似乎是问题的根源——如果你不知道并行指令的存在，你就会就像提问者对“更多”如何神奇地变成“更少”一样感到困惑。不过，依赖性和瓶颈——编译器如何确定哪些优化选项可用——将是很好的后续问题。 (24认同)
- 但指令级并行性与 SIMD 并行性同样重要。也许更重要的是，每个向量只有 2 个“double”，而 SIMD FP“addsd”/“addpd”在 Nehalem 和 Sandy Bridge 上具有 3 个周期延迟、1 个周期吞吐量。（尽管循环中有两个独立的加法链，这可能会导致每 1.5 个时钟周期进行一个标量 FP 加法，所以是的，SIMD 可能更重要。） (13认同)
- 无论如何，跨循环迭代的**串行依赖**实际上*是*并行与串行代码（以及该代码的执行）的最终关键，并且IMO将是一个很好的开头段落。编译器和 CPU 可以通过多种方式利用它，编译器自动向量化，CPU 利用独立循环迭代的 ILP。即使您只想讨论 SIMD 矢量化，发现循环中可用的数据并行性也是关键的第一个观察结果。（我确实已经赞成这个答案；总体不错，但如果它是从并行性与依赖关系开始的话，我会更喜欢它） (12认同)
- 在您的更新中，您提到了[强度减少优化](https://en.wikipedia.org/wiki/Strength_reduction)。问题中提出的优化*是*一种强度降低的奇特案例，用循环进位加法链代替独立乘法。因此，如果编译器这样做（使用“-ffast-math”），您希望它以一种展开友好的方式进行，以允许矢量化。 (8认同)
- @Nimitz14不过，您实际上确实指出了一个好点：这些是双精度数，而不是浮点数，因此每个寄存器仅保存_two_，而不是_four_，并且每个循环迭代的总数是_four_，而不是_8_。我的错误，我相应地更新了答案。接得好！ (2认同)
Chr*_*odd 65
这里的主要区别是循环依赖性。第二种情况下的循环依赖于\xe2\x80\x94 循环中的操作依赖于先前的迭代。这意味着每次迭代在上一次迭代完成之前都无法开始。在第一种情况下，循环体是完全独立的\xe2\x80\x94循环体中的所有内容都是自包含的，仅取决于迭代计数器和常量值。这意味着循环可以并行计算\xe2\x80\x94多个迭代可以同时工作。然后，这允许循环被简单地展开和矢量化，从而重叠许多指令。
\n如果您要查看性能计数器（例如，使用perf stat ./1），您会发现第一个循环除了运行速度更快之外，每个周期还运行更多指令 (IPC)。相比之下，第二个循环有更多的依赖周期\xe2\x80\x94时间，此时CPU什么都不做，等待指令完成，然后才能发出更多指令。
第一个可能会成为内存带宽瓶颈，特别是如果您让编译器在Sandy Bridge ( gcc -O3 -march=native) 上使用 AVX 自动矢量化（如果它设法使用 256 位矢量）。那时，IPC 将下降，特别是对于L3 缓存来说输出数组太大。
\n
需要注意的是：展开和矢量化不需要独立的循环\xe2\x80\x94，当（某些）循环依赖项存在时，您可以执行它们。然而，它更难，回报也更少。因此，如果您想看到矢量化的最大加速，它有助于尽可能消除循环依赖性。
\n- 为了通过手动向量化获得*最大*收益，我认为您实际上可能会使用版本 2，但是使用多个以锁步方式前进的向量，四个不同的 Z 和 Y 向量，例如“Z0 += 8*A2”（或如果每个 Z 向量包含 4 个双精度数而不是 2 个，则为“16*A2”）。您需要一些数学来解释将元素跨过 8 或 16 个“i”值而不是 1，也许是在某个地方进行乘法。几乎可以肯定，您可以比每次迭代重做 int->FP 转换做得更好；这是获得独立迭代的昂贵方法。 (8认同)
Pet*_*des 42
这种有限差分强度降低优化方法可以比单独重新评估每个 多项式的最佳方法提供加速i。但只有当你将其推广到更大的步幅时，循环中仍然具有足够的并行性。我的版本在Skylake上每个时钟周期存储一个向量（四个双精度），用于适合 L1d 缓存的小型数组；否则就是带宽测试。在早期的 Intel 上，它还应该最大化 SIMD FP 添加吞吐量，包括带有 AVX 的Sandy Bridge（如果对齐输出，则每两个时钟 1 个 256 位添加/时钟，以及每两个时钟 1 个 256 位存储。）
对上一次迭代的值的依赖是致命的
这种强度降低优化（只是相加而不是从新的i和相乘开始）引入了跨循环迭代的串行依赖性，涉及 FP 数学而不是整数增量。
原始版本在每个输出元素上都具有数据并行性：每个输出元素仅取决于常量及其自身的i值。编译器可以使用 SIMD（ SSE2或AVX ，如果使用的话）自动矢量化-O3 -march=native，CPU 可以通过乱序执行来重叠循环迭代的工作。尽管有大量额外工作，CPU 仍能够在编译器的帮助下应用足够的强力。
poly(i+1)但以计算的版本poly(i)并行性非常有限；没有 SIMD 矢量化，并且您的 CPU 每四个周期只能运行两个标量加法，例如，其中四个周期是 Intel Skylake 到 Tiger Lake 上 FP 加法的延迟。（https://uops.info/）。
huseyin tugrul buyukisik 的答案展示了如何在更现代的 CPU 上接近最大化原始版本的吞吐量，使用两个 FMA 操作来评估多项式（霍纳方案），再加上整数到浮点转换或浮点数点增量。（后者创建了一个 FP 依赖链，您需要展开才能隐藏它。）
因此，最好的情况是每个SIMD输出向量有 3 个浮点数学运算。（加上一家商店）。当前的 Intel CPU 仅具有两个浮点执行单元，可以运行包括 int 到 double 在内的 FP 数学运算。（对于 512 位向量，当前的 CPU 会关闭端口 1 上的向量 ALU，因此只有两个 SIMD ALU 端口，因此非 FP 数学运算（例如 SIMD 整数增量）也将争夺 SIMD 吞吐量。除了对于只有一个 512 位 FMA 单元的 CPU，则端口 5 可用于其他工作。）
自 Zen 2 起，AMD 在两个端口上有两个 FMA/mul 单元，在两个不同端口上有两个 FP add/sub 单元，因此，如果您使用 FMA 进行加法，理论上每个时钟周期最多有四个 SIMD 加法。
Haswell/Broadwell 有 2/时钟 FMA，但只有 1/时钟 FP 加/减（延迟较低）。这对于简单的代码来说是有好处的，但对于经过优化以具有大量并行性的代码来说则不太好。这可能就是英特尔在 Skylake 中改变它的原因。（Alder Lake 重新引入了低延迟 FP 加/减，但 2/时钟吞吐量与乘法相同。有趣的是，非交换延迟：目的地只有 2 个周期，其他操作数有 3 个周期，因此对于较长的依赖链来说非常有用.)
您的 Sandy Bridge (E5-1620) 和 Nehalem (W5580) CPU 在单独的端口上具有 1/时钟 FP 加/减、1/时钟 FP 乘。这就是 Haswell 的基础。为什么添加额外的乘法不是一个大问题：它们可以与现有的添加并行运行。（Sandy Bridge 的宽度为 256 位，但您在编译时未启用 AVX：使用-march=native。）
寻找并行性：以任意步幅减少强度
您compute2根据前一个值计算下一个 Y 和下一个 Z。即，步幅为 1 时，您需要的值data[i+1]。因此，每次迭代都依赖于前一次迭代。
如果将其推广到其他步幅，则可以将 Y 和 Z 值推进 4、6、8 或更多单独的 Y 和 Z 值，以便它们彼此同步跳跃，并且彼此独立。这将为编译器和/或 CPU 重新获得足够的并行性以供利用。
poly(i) = A i^2 + B i + C
poly(i+s) = A (i+s)^2 + B (i+s) + C
= A*i^2 + A*2*s*i + A*s^2 + B*i + B*s + C
= poly(i) + A*2*s*i + A*s^2 + B*s + C
所以这有点混乱，不完全清楚如何将其分解为 Y 和 Z 部分。（这个答案的早期版本弄错了。）
可能更容易从一阶和二阶差分向后推算 FP 值序列的跨步（有限差分法）。这将直接找到我们需要添加的内容；Z[] 初始化器和步骤。
这基本上就像取一阶和二阶导数，然后优化循环有效积分以恢复原始函数。以下输出是由main下面基准测试的正确性检查部分生成的。
# method of differences for stride=1, A=1, B=0, C=0
poly(i) 1st 2nd difference from this poly(i) to poly(i+1)
0 1
1 3 2 # 4-1 = 3 | 3-1 = 2
4 5 2 # 9-4 = 5 | 5-3 = 2
9 7 2 # ...
16 9 2
25 11 2
相同的多项式 ( x^2)，但以 3 的步幅求差。非 2 的幂有助于显示步幅的因子/幂来自何处，与自然发生的因子 2 相比。
# for stride of 3, printing in groups. A=1, B=0, C=0
poly(i) 1st 2nd difference from this poly(i) to poly(i+3)
0 9
1 15
4 21
9 27 18 # 36- 9 = 27 | 27-9 = 18
16 33 18 # 49-16 = 33 | 33-15 = 18
25 39 18 # ...
36 45 18 # 81-36 = 45 | 45-27 = 18
49 51 18
64 57 18
81 63 18
100 69 18
121 75 18
Y[] 和 Z[] 初始值设定项
初始值
Y[j] = poly(j)，因为它必须存储到相应位置 (data[i+j] = Y[j]) 的输出中。初始的
Z[j]将被添加到Y[j]，并且需要将其添加到poly(j+stride). 因此，初始Z[j] = poly(j+stride) - Y[j]，如果需要，我们可以对其进行代数简化。（对于编译时常量 A、B、C，编译器将以任一方式传播常量。）Z[j]保持跨越 的一阶差异poly(x)，作为 的起点poly(0..stride-1)。这是上表中的中间列。必要的更新
Z[j] += second_difference是一个标量常数，正如我们从二阶差分中看到的那样，它是相同的。通过使用几个不同的
stride和值（ i^ 2A的系数），我们可以看到它是。（使用 3 和 5 等非互质值有助于理清问题。）通过更多的代数，您可以象征性地展示这一点。从微积分 PoV 来看，因子 2 有意义：，二阶导数为。A * 2 * (stride * stride)d(A*x^2)/dx = 2Ax2A
poly(i) = A i^2 + B i + C
poly(i+s) = A (i+s)^2 + B (i+s) + C
= A*i^2 + A*2*s*i + A*s^2 + B*i + B*s + C
= poly(i) + A*2*s*i + A*s^2 + B*s + C
对于stride=1（不展开）这些简化为原始值。但对于较大的stride，编译器可以将 Y[] 和 Z[] 的元素分别保留在几个 SIMD 向量中，因为每个向量Y[j]仅与相应的 进行交互Z[j]。
编译器 (SIMD) 和 CPU（流水线执行单元）可以stride利用独立的并行依赖链，运行stride时间比原始版本更快，直到您在 SIMD FP 添加吞吐量（而不是延迟）上遇到瓶颈，或者存储带宽（如果您的缓冲区不适合 L1d）。（或者直到编译器面对这些循环并且不很好地/根本不展开和矢量化这些循环！）
在实践中如何编译：与 clang 配合得很好
（Godbolt 编译器资源管理器）使用 clang14 很好地进行 Clang 自动向量化stride=16（4x YMM 向量，double每个向量 4 秒）-O3 -march=skylake -ffast-math。
看起来 clang 进一步展开了 2，快捷方式Z[j] += diff2进入tmp = Z[j] + diff2;/ Z[j] += 2*diff2;。这减轻了 Z 依赖链的压力，只留下 Y[j] 来应对 Skylake 上的延迟瓶颈。
因此，每个 asm 循环迭代都会执行 2x 8vaddpd条指令和 2x 4 条存储。循环开销是add+ 宏融合cmp/jne，所以 2 uop。（或者对于全局数组，只有一个add/jne微指令，将负索引向上计数到零；它相对于数组末尾进行索引。）
vaddpdSkylake 在每个时钟周期几乎 1 次存储和 2 次运行。这是这两件事的最大吞吐量。前端只需要跟上略高于 3 uops/时钟周期的速度，但自 Core2 以来已经是 4 宽了。Sandy Bridge 系列中的 uop 缓存使这不再是问题。（除非你在 Skylake 上遇到 JCC 勘误表，所以我曾经-mbranches-within-32B-boundaries 使用 clang pad 指令来避免这种情况。）
由于 Skylake 的vaddpd延迟为 4 个周期，4 个依赖链stride=16勉强足以保持 4 个独立操作的运行。每当 aY[j]+=不运行它准备好的循环时，就会产生泡沫。由于 Clang 对 Z[] 链的额外展开，Z[j]+= 可以提前运行，因此 Z 链可以领先。通过最旧就绪优先调度，它往往会稳定在 Yj+= uop 不存在冲突的状态，显然，因为它确实在我的 Skylake 上全速运行。如果我们能让编译器仍然为 生成良好的 asm stride=32，那将留下更多空间，但不幸的是它没有。（以奇怪尺寸的更多清理工作为代价。）
奇怪的是，Clang 只用 向量化它-ffast-math。下面完整基准测试中的模板版本不需要-ffast-math. 源代码经过精心编写，适合 SIMD，并按源代码顺序进行数学运算。（不过，快速数学使得 clang 能够更多地展开 Z 增量。）
编写循环的另一种方法是使用一个内部循环而不是所有 Y 操作，然后是所有 Z 操作。这在下面的基准测试中很好（有时实际上更好），但这里即使使用-ffast-math. 对于像这样的重要问题，从编译器中获得最佳的展开 SIMD 汇编可能会很繁琐且不可靠，并且可能需要一些尝试。
我将它包含在Godbolt 上的#if 0//块中。#else#endif
# method of differences for stride=1, A=1, B=0, C=0
poly(i) 1st 2nd difference from this poly(i) to poly(i+1)
0 1
1 3 2 # 4-1 = 3 | 3-1 = 2
4 5 2 # 9-4 = 5 | 5-3 = 2
9 7 2 # ...
16 9 2
25 11 2
我们必须手动选择合适的展开金额。太大的展开因子甚至会阻止编译器了解正在发生的情况并停止将临时数组保留在寄存器中。例如，32or24是 clang 的问题，但不是16。可能有一些调整选项可以强制编译器将循环展开到一定数量；GCC 有时候可以用来让它在编译时看透一些东西。
另一种方法是使用#include <immintrin.h>和__m256d Z[4]代替进行手动矢量化double Z[16]。但这个版本可以针对其他 ISA（例如 AArch64）进行矢量化。
大展开因子的其他缺点是，当问题大小不是展开因子的倍数时，会留下更多的清理工作。（您可以使用该compute1策略进行清理，让编译器在执行标量之前对其进行矢量化一两次迭代。）
理论上，编译器可以通过对原始多项式进行强度降低或查看步幅如何累积来为您执行此操作。-ffast-mathcompute1compute2
但在实践中，这确实很复杂，而且人类必须自己做。除非/直到有人开始教编译器如何寻找这样的模式并应用差异方法本身，并选择步幅！但是，即使使用-ffast-math. （整数不会有任何精度问题，但它仍然是一个复杂的模式匹配/替换。）
实验性能结果：
我在我的台式机（i7-6700k）上使用 clang13.0.0 进行了测试。事实上，它确实在每个时钟周期运行 1 个 SIMD 存储，并具有多种编译器选项组合（无论是否是快速数学）以及#if 0内部#if 1循环策略。我的基准/测试框架基于 @huseyin tugrul buyukisik 的版本，经过改进以在rdtsc指令之间重复更多可测量的量，并使用测试循环来检查多项式的简单计算的正确性。
我还让它补偿了核心时钟频率和TSC 读取的“参考”频率rdtsc之间的差异，在我的例子中是 3.9 GHz 与 4008 MHz。（额定最大睿频为 4.2 GHz，但balance_performance在 Linux 上使用 EPP = 时，它只想将时钟频率提高到 3.9 GHz。）
Godbolt 上的源代码：使用一个内部循环，而不是 3 个单独的j<16循环，并且不使用-ffast-math. 用于__attribute__((noinline))防止其内联到重复循环中。选项和源的其他一些变化导致了vpermpd循环内的一些混乱。
下面的基准数据来自以前的版本，带有错误的 Z[j] 初始值设定项，但循环 asm 相同。Godbolt 链接现在在定时循环之后进行了正确性测试，该测试通过了。double即使没有#if 1/-ffast-math来允许 clang 额外展开，实际性能在我的桌面上仍然相同，每个 只超过 0.25 个周期。
# for stride of 3, printing in groups. A=1, B=0, C=0
poly(i) 1st 2nd difference from this poly(i) to poly(i+3)
0 9
1 15
4 21
9 27 18 # 36- 9 = 27 | 27-9 = 18
16 33 18 # 49-16 = 33 | 33-15 = 18
25 39 18 # ...
36 45 18 # 81-36 = 45 | 45-27 = 18
49 51 18
64 57 18
81 63 18
100 69 18
121 75 18
fp_arith_inst_retired.256b_packed_double对于每个 FP add 或 mul 指令计数为 1（对于 FMA 为 2），因此整个程序的每个时钟周期我们获得 1.98vaddpd条指令，包括打印等。这非常接近理论最大 2/时钟，显然没有受到次优微指令调度的影响。（我增加了重复循环，因此程序将大部分总时间花在此处，从而使perf stat整个程序变得有用。）
此优化的目标是以更少的 FLOPS 完成相同的工作，但这也意味着我们在不使用 FMA 的情况下实质上达到了 Skylake 8 FLOP/时钟的限制。（单核上 3.9 GHz 时为 30.58 GFLOP/s）。
非内联函数的汇编 ( objdump -drwC -Mintel); clang 使用 4 个 Y、Z 对的 YMM 向量，并将循环进一步展开 3 倍，使其成为 24 KiB 大小的精确倍数，且无需清理。请注意add rax,0x30每次迭代都会执行 3 * stride=0x10 双倍操作。
0000000000001440 <void compute2<3072>(double*)>:
# just loading constants; the setup loop did fully unroll and disappear
1440: c5 fd 28 0d 18 0c 00 00 vmovapd ymm1,YMMWORD PTR [rip+0xc18] # 2060 <_IO_stdin_used+0x60>
1448: c5 fd 28 15 30 0c 00 00 vmovapd ymm2,YMMWORD PTR [rip+0xc30] # 2080
1450: c5 fd 28 1d 48 0c 00 00 vmovapd ymm3,YMMWORD PTR [rip+0xc48] # 20a0
1458: c4 e2 7d 19 25 bf 0b 00 00 vbroadcastsd ymm4,QWORD PTR [rip+0xbbf] # 2020
1461: c5 fd 28 2d 57 0c 00 00 vmovapd ymm5,YMMWORD PTR [rip+0xc57] # 20c0
1469: 48 c7 c0 d0 ff ff ff mov rax,0xffffffffffffffd0
1470: c4 e2 7d 19 05 af 0b 00 00 vbroadcastsd ymm0,QWORD PTR [rip+0xbaf] # 2028
1479: c5 fd 28 f4 vmovapd ymm6,ymm4 # buggy Z[j] initialization in this ver used the same value everywhere
147d: c5 fd 28 fc vmovapd ymm7,ymm4
1481: c5 7d 28 c4 vmovapd ymm8,ymm4
1485: 66 66 2e 0f 1f 84 00 00 00 00 00 data16 cs nop WORD PTR [rax+rax*1+0x0]
# top of outer loop. The
-
@huseyintugrulbuyukisik：使用我的 Skylake 桌面的测试结果进行更新：它确实每个时钟运行 1 个存储（以及两个“vaddpd”），因此在没有 AVX-512（我的桌面没有）的情况下，每个元素得到 0.251 个周期。在测试时，我注意到您使用的是“rdtsc”数字而不是核心时钟周期，这是一个很大的假设。对于某些 Xeon 来说，当[运行“重”512 位指令](/sf/ask/3979696871/) 时，实际核心时钟接近 TSC 频率可能是正确的。频率），但这是一个有风险的假设。 (2认同)
-
但无论如何，大概与我的 asm 相同但使用 ZMM 向量也可以在 Skylake-avx512 CPU 上每个时钟运行 1 个存储，因此每个元素大约 0.125 个周期。如果没有覆盖调整启发式的选项，让编译器生成这样的 asm 可能会出现问题，因此如果您不使用内在函数，则可能会出现潜在的实际问题。 (2认同)
-
@huseyintugrulbuyukisik：我们并不知道代码最终运行的服务器实例的 CPU 频率，尽管我们可以使用 CPUID 来获取品牌字符串并打印它，其中可能包括库存“额定”频率。有了它，就可以进行手动计算（或修正 RDTSC 猜测数字）。或许可以采用 Quick-bench 的 NOP 循环计时策略来估计当前 CPU 频率，尽管运行 AVX-512“重”指令导致的睿频降低使这变得更加困难。 (2认同)
-
无论如何，这只是一个理论问题；没有必要为生产用途而实际优化它，只是概念证明就可以了。因此，从纯 C++ 源代码自动矢量化并不是我要花更多时间的事情，除非/除非在特定项目中出现实际用例，该项目将控制我们使用哪些编译器/选项。可以使用什么，要调整什么问题大小，以及如何调用等等。 (2认同)
-
@huseyintugrulbuyukisik：是的，即使对于以前版本的算法，在许多情况下也是如此。除非您想在 ALU 吞吐量瓶颈的循环中多次重新读取它，否则可能值得保留。（特别是如果您可以缓存块它，这样您就不会浪费系统范围的内存带宽，或者如果您的其他循环也需要的话，则可以浪费 L3 或 L2 带宽。） (2认同)
gna*_*729 12
如果您需要此代码快速运行，或者您很好奇，您可以尝试以下操作：
您将 a[i] = f(i) 的计算更改为两次加法。修改它以使用两次加法计算 a[4i] = f(4i)，使用两次加法计算 a[4i+1] = f(4i+1)，依此类推。现在您有四个可以并行完成的计算。
编译器很有可能会执行相同的循环展开和向量化，并且您有相同的延迟，但对于四个操作，而不是一个。
hus*_*sik 10
通过仅使用加法作为优化，您会丢失（较新 CPU 的）乘法管道的所有GFLOPS，并且循环携带的依赖性会通过停止自动矢量化而使情况变得更糟。如果它是自动向量化的，它会比乘法+加法快得多。每个数据的能源效率更高（仅相加比相加更好）。
另一个问题是，由于累积的加法次数，数组末尾会收到更多舍入错误。但它不应该在非常大的数组之前可见（除非数据类型变为浮点型）。
当您应用带有 GCC 构建选项的Horner 方案时（在较新的 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 优化手册中的优化一定来自 Quake 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 数组按数据数组大小对齐且足够小时，它以每个元素0.26 个周期alignas(64)的速度运行。
- 询问者仅在 Nehalem 和 Sandybridge Xeon CPU 上进行测试，这些 CPU 不支持 FMA。您忘记提及用于让它使用 AVX2+FMA 自动矢量化的构建选项。但是，是的，如果您有 FMA，这是一个很好的策略。即使您不这样做，在 CPU 上，“mulpd”也可能在与“addpd”不同的端口上运行，因此它们只竞争前端吞吐量。如果你只关心速度，而不是准确性，那么 gnasher 的答案中建议的策略（我之前在评论中建议过）使用多个累加器来隐藏 FP 延迟，可能仍然更好，避免了 int->FP 成本。 (4认同)
与加法相比，乘法在我们的 CPU 中以显着慢着称。
这在历史上可能是正确的，对于更简单的低功耗 CPU 来说可能仍然如此，但如果 CPU 设计者准备好“解决问题”，乘法几乎可以像加法一样快地进行。
现代 CPU 被设计为通过流水线和多个执行单元的组合来同时处理多个指令。
但这样做的问题是数据依赖性。如果一条指令依赖于另一条指令的结果，那么在它所依赖的指令完成之前，它的执行无法开始。
现代 CPU 试图通过“乱序执行”来解决这个问题。等待数据的指令可以保持排队，同时允许执行其他指令。
但即使采取了这些措施，有时 CPU 也可能会耗尽要调度的新工作。
- 对于从 Skylake 开始（Alder Lake 之前）的 Intel CPU 上的 FP 来说，情况确实如此。FP add/sub/mul/fma 都具有完全相同的性能，在相同的 2 个（完全流水线）执行端口上运行，具有相同的 4 周期延迟。Alder Lake 将 FP add/sub 加速到 3 个周期，就像在 Haswell 中一样（但仍然具有像 mul/fma 一样的 2/时钟吞吐量，[与 Haswell 不同](https:// electronics.stackexchange.com/questions/452181/why -intels-haswell-chip 是否允许浮点乘法两次/452366#452366))。 (6认同)
- 但对于整数数学则不然。对于标量整数，具有 3 个周期延迟的 1/时钟与具有 1c 延迟的 4/时钟相比，现代 Intel 上的 SIMD 整数吞吐量也是 4 倍。与旧 CPU 相比，整数乘法的吞吐量仍然相当高。 (3认同)
| 归档时间： |
|
| 查看次数： |
86631 次 |
| 最近记录： |