我正在閱讀 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 次加法！
所以我繼續添加一個回圈來執行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）。
我承認我沒想到會有這樣的結果。不是因為我是 Agner 的粉絲（我是，但這無關緊要）。
知道我缺少什么嗎？我在這里犯了什么錯誤，可以解釋速度的差異嗎？請注意，我已經在 Xeon W5580 和 Xeon E5 1620 上進行了測驗——在這兩個版本中，第一個（啞）版本比第二個版本快得多。
編輯：為了簡單地再現結果，我在兩個版本的代碼中添加了兩個要點：Dumb 但不知何故更快和優化，但不知何故更慢。
PS請不要評論浮點精度問題；這不是這個問題的重點。
uj5u.com熱心網友回復：
了解您所看到的性能差異的關鍵在于矢量化。是的，基于加法的解決方案在其內部回圈中只有兩條指令，但重要的區別不在于回圈中有多少指令，而在于每條指令執行了多少作業。
在第一個版本中，輸出完全依賴于輸入：每個函式本身data[i]就是一個函式i，這意味著每個函式都data[i]可以按任何順序計算：編譯器可以向前、向后、橫向等等，你會仍然得到相同的結果 - 除非您從另一個執行緒觀察該記憶體，否則您永遠不會注意到資料正在以哪種方式被處理。
在第二個版本中，輸出不依賴于i——它依賴于回圈的最后一次的A和。Z
如果我們將這些回圈的主體表示為小的數學函式，它們將具有非常不同的整體形式：
- f(i) -> 迪
- f(Y, Z) -> (di, Y', Z')
在后一種形式中，沒有實際的依賴關系i——計算函式值的唯一方法是知道函式的前一次Y和Z最后一次呼叫，這意味著函式形成一個鏈——你不能這樣做下一個，直到你完成上一個。
為什么這很重要？因為 CPU 有向量并行指令，每條指令可以同時執行四個算術運算！（AVX CPU 可以并行執行更多操作。）即四次乘法、四次加法、四次減法、四次比較——四次任意！因此，如果您嘗試計算的輸出僅取決于輸入，那么您可以安全地一次執行四個 - 無論它們是向前還是向后，因為結果是相同的。但是，如果輸出依賴于先前的計算，那么您將被困在串行形式中——一次一個。
這就是為什么“更長”的代碼會贏得性能。盡管它有更多的設定，并且實際上做了更多的作業，但大部分作業都是并行完成的：它不僅僅是data[i]在回圈的每次迭代中計算——它是計算data[i], data[i 1], data[i 2],data[i 3]同時，然后跳到下一組四個。
為了擴展我在這里的意思，編譯器首先將原始代碼變成了這樣的東西：
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 有特殊的內置指令，可以對不同的資料執行四*到四個 操作同時，每個僅在一個時鐘周期內。
請注意p較快解決方案中說明中的字母 -addpd和mulpd- 以及s較慢解決方案中說明中的字母 - addsd。那是“添加包裝雙打”和“乘以包裝雙打”與“添加單雙打”。
不僅如此，看起來編譯器也部分展開了回圈——回圈不僅每次迭代執行四個值，實際上是八個，并且交錯操作以避免依賴和停頓，這減少了回圈的次數匯編代碼也必須進行測驗i < 1000。
但是，所有這些都只有在回圈的迭代之間沒有依賴關系的情況下才有效：如果唯一決定每個回圈發生什么的事情data[i]就是i它自己。如果存在依賴關系，如果上一次迭代的資料影響下一次迭代，那么編譯器可能會受到它們的限制，以至于它根本無法更改代碼——而不是編譯器能夠使用花哨的并行指令或聰明的優化（CSE、強度降低、回圈展開、重新排序等），你得到的代碼正是你輸入的代碼——添加 Y，然后添加 Z，然后重復。
但是在這里，在代碼的第一個版本中，編譯器正確地識別出資料中沒有依賴關系，并發現它可以并行完成作業，所以它確實做到了，這就是一切的不同之處。
uj5u.com熱心網友回復：
這里的主要區別是回圈依賴。第二種情況下的回圈是依賴的——回圈中的操作依賴于之前的迭代。這意味著每次迭代甚至不能在前一次迭代完成之前開始。在第一種情況下，回圈體是完全獨立的——回圈體中的所有內容都是自包含的，僅取決于迭代計數器和常量值。這意味著回圈可以并行計算——多個迭代可以同時作業。然后，這允許回圈被簡單地展開和矢量化，重疊許多指令。
如果您要查看性能計數器（例如，使用perf stat ./1），您會看到第一個回圈除了運行速度更快之外，每個回圈還運行更多指令（IPC）。相比之下，第二個回圈有更多的依賴周期——CPU 無所事事，等待指令完成，然后才能發出更多指令的時間。
第一個可能是記憶體帶寬的瓶頸，特別是如果您讓編譯器在 Sandybridge ( gcc -O3 -march=native) 上使用 AVX 自動矢量化，如果它設法使用 256 位向量。那時 IPC 將下降，特別是對于 L3 快取來說太大的輸出陣列。
請注意，展開和矢量化不需要獨立的回圈——當（某些）回圈依賴項存在時，您可以執行它們。然而，它更難，回報也更少。因此，如果您想從矢量化中看到最大的加速，它有助于在可能的情況下消除回圈依賴。
uj5u.com熱心網友回復：
如果您需要此代碼快速運行，或者您對此感到好奇，可以嘗試以下方法：
您將 a[i] = f(i) 的計算更改為兩個加法。修改它以使用兩次加法計算 a[4i] = f(4i)，使用兩次加法計算 a[4i 1] = f(4i 1)，依此類推。現在您有四個可以并行完成的計算。
編譯器很有可能會執行相同的回圈展開和矢量化，并且您具有相同的延遲，但是對于四個操作，而不是一個。
uj5u.com熱心網友回復：
通過僅使用加法作為優化，您將丟失（較新 CPU 的）乘法管道的所有 gflops，并且回圈攜帶的依賴性通過停止自動矢量化使其變得更糟。
另一個問題是由于累加的數量，陣列的末尾會收到更多的舍入誤差。但它不應該在非常大的陣列之前可見（除非資料型別變為浮點數）。
當您將 Horner Scheme 與 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 中，所有 FP 操作都是“打包”向量形式和更少的指令（它是兩次展開的版本）。每 64 位元組資料有 16 條指令（如果是 AVX512，則為 128 位元組）。
Horner Scheme 的另一個好處是它在 FMA 指令中的計算精度更高，并且每次回圈迭代只有 O(1) 次操作，因此它不會在更長的陣列中累積太多錯誤。
我認為 Agner Fog 優化手冊中的優化必須來自 Quake-3 快速反平方根近似時代。在那個時候，SIMD 不夠寬，無法產生太大的影響，也缺乏對 sqrt 函式的支持。手冊說著作權所有 2004，所以賽揚有 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};
// 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
沒有不可隱藏的 int-to-double 轉換延遲，沒有回圈依賴，精度更高，每 128 位元組資料有 13 條指令（每個double型別資料元素約 0.8 條指令）
當資料、A、B 和 C 陣列按足夠小的資料陣列大小對齊時alignas(64)，它以每個元素速度0.26 個周期運行。
uj5u.com熱心網友回復：
與加法相比，乘法在我們的 CPU 中速度明顯較慢而聞名。
這在歷史上可能是正確的，并且對于更簡單的低功耗 CPU 可能仍然正確，但是如果 CPU 設計人員準備“解決問題”，則乘法幾乎可以與加法一樣快。
現代 CPU 旨在通過流水線和具有多個執行單元的組合同時處理多條指令。
但是，問題在于資料依賴性。如果一條指令依賴于另一條指令的結果，則在它所依賴的指令完成之前無法開始執行。
現代 CPU 試圖通過“亂序執行”來解決這個問題。等待資料的指令可以保持排隊，同時允許執行其他指令。
但即使采取了這些措施，有時 CPU 也會簡單地用完新的作業來調度。
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/479064.html