1.3.2　積分計算圓周率π

這個簡單的例子展示了使用指令級并行、數據級并行技術的初級方法，并且不涉及存儲器和緩存的使用。

使用以下公式：

并采用離散積分的方法求π值。

1）給出一個簡單的版本1作為基準（baseline），具體如代碼清單1-7所示。

代碼清單1-7　版本1：計算圓周率的基準版

---

double compute_pi_naive(size_t dt)  // version 1
{
   double pi = 0.0;
   double delta = 1.0 / dt;
   for (size_t i = 0; i < dt; i++) {
   　double x = (double)i / dt;
   　pi += delta / (1.0 + x * x);
    }
    return pi * 4.0;
}

---

dt為[0，1]區間分段的大小，理論上dt越大計算結果越精確，但同時計算量也越大；delta為每個分段的長度；x是第i個分段的橫坐標，用它帶入積分函數，并乘以delta表示每個分段的面積；用pi累加所有的分段面積，乘以4并返回結果。

在筆者的Intel i3 2100 SNB處理器上，這個版本在輸入dt=128M時的執行時間約為1.919秒。為了方便比較，將這個版本的性能標準化，記為1。

2）下面的版本2將在單線程內使用SIMD指令數據并行技術進一步縮短執行時間。SNB處理器支持AVX指令集，這里使用AVX指令來向量化單線程內部的計算，如代碼清單1-8所示。

代碼清單1-8　版本2：在線程內部使用SIMD數據并行技術加速線程執行

---

double compute_pi_omp_avx(size_t dt){
    double pi = 0.0;
    double delta = 1.0 / dt;
    __m256d ymm0, ymm1, ymm2, ymm3, ymm4;
    ymm0 = _mm256_set1_pd(1.0);
    ymm1 = _mm256_set1_pd(delta);
    ymm2 = _mm256_set_pd(delta * 3, delta * 2, delta, 0.0);
    ymm4 = _mm256_setzero_pd();
    for (int i = 0; i <= dt-4; i += 4) {
   　　 ymm3 = _mm256_set1_pd(i * delta);
   　　 ymm3 = _mm256_add_pd(ymm3, ymm2);
   　　 ymm3 = _mm256_mul_pd(ymm3, ymm3);
   　　 ymm3 = _mm256_add_pd(ymm0, ymm3);
   　　 ymm3 = _mm256_div_pd(ymm1, ymm3);
   　　 ymm4 = _mm256_add_pd(ymm4, ymm3);
    }
    double tmp[4] __attribute__((aligned(32)));
    _mm256_store_pd(tmp, ymm4);
    pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
    return pi * 4.0;
}

---

_mm256_set1_pd（）語句用于將一個double型浮點值賦給一個_m256變量的所有位置；_mm256_set_pd（e3，e2，e1，e0）是將4個參數依次放入_m256變量的各個位置，參數順序和存儲順序相反。注意，這兩個intrinsic實際上并不對應單一一條指令，而是多條指令混合而來，故應盡量減少將其用于最內層循環。

程序的最內層循環的前兩條語句首先構造出遞增的x的向量表示，然后進行向量的自乘加1，再被delta的向量除，加到一個累加向量ymm4中。

這樣就完成了初步的向量化工作。經測試，同樣在i32100的CPU上，相比沒做向量化的程序加速比達到2.6x左右。256位的向量寄存器可以同時存儲計算4個double類型的浮點數，但由于浮點標量除法和向量除法的延遲不同，并不能達到4x的加速比。

3）初步的向量化并沒有充分填充向量指令的流水線，即沒有掩蓋浮點計算的延遲。版本3通過使用循環展開技術，以使用更多的寄存器，來解除數據依賴，填充浮點計算流水線（這里僅列出循環內代碼），如代碼清單1-9所示。

代碼清單1-9　版本3：使用循環展開技術解除數據依賴

---

ymm5 = _mm256_setzero_pd();
for (i = 0; i <= dt - 8; i += 8){
   　　ymm3 = _mm256_set1_pd(i * delta);
   　　ymm3 = _mm256_add_pd(ymm3, ymm2);
   　　ymm3 = _mm256_mul_pd(ymm3, ymm3);
   　　ymm3 = _mm256_add_pd(ymm0, ymm3);
   　　ymm3 = _mm256_div_pd(ymm1, ymm3);
   　　ymm4 = _mm256_add_pd(ymm4, ymm3);
   　　ymm6 = _mm256_set1_pd((i+4)* delta);
   　　ymm6 = _mm256_add_pd(ymm6, ymm2);
   　　ymm6 = _mm256_mul_pd(ymm3, ymm6);
   　　ymm6 = _mm256_add_pd(ymm0, ymm6);
   　　ymm6 = _mm256_div_pd(ymm1, ymm6);
   　　ymm5 = _mm256_add_pd(ymm5, ymm6); 
}
ymm4 = _mm256_add_pd(ymm4, ymm5);

---

ymm3和ymm4串起兩條獨立的指令依賴關系，在X86這種亂序多發射的指令調度里可以并行執行，或者互相填充流水線，最后都加到ymm5上。測試表明，相對初步向量化加速50%左右。

4）將前面3個版本綜合，如表1-9所示。

表1-9　積分求π（3個版本性能比較）

從表中可以看出，SIMD指令向量化和循環展開以優化流水線效率，都獲得了比較好的性能提升。由于本章的主要目的是介紹Intel X86 CPU上的SSE/AVX指令優化技術，故并不會提到多線程優化技術。