如何優化這個數學運算的速度

Question

我試圖優化一個功能，佔用大量的執行時間，它會多次計算下面的數學運算。無論如何要讓這個操作更快？

float total = (sqrt(
      ((point_A[j].length)*(point_A[j].length))+ 
      ((point_B[j].width)*(point_B[j].width))+ 
      ((point_C[j].height)*(point_C[j].height)) 
                 )); 

Answer 1

通常，您希望避免使用傳統的幾何圖形和三角函數，只要有意義就切換到矢量運算。例如這意味着以平方長度而不是長度來工作。許多使用長度的算法可以很容易地修改，以便使用平方長度。但是，如果您必須採取平方根，我會建議在您的情況下嘗試使用sqrt(x*x + y*y)的專用函數hypot(x,y)（在這裏，您必須調用它兩次：例如hypot(x,hypot(y,z))）。這可能或不會幫助。

另外，還要考慮sqrtf代替sqrt，並接通編譯器優化更快數學（例如-ffast-math爲gcc）或優化（或不同的庫），該犧牲精度速度。

Answer 2

如果內存很便宜，那麼您可以執行以下操作，從而提高命中率CPU cache。既然你沒有發佈更多的細節，所以我會在這裏做一些假設。

long tmp_len_square[N*3]; 

for (int j = 0; j < N; ++j) { 
    tmp_len_square[3 * j] = (point_A[j].length)*(point_A[j].length); 
} 

for (int j = 0; j < N; ++j) { 
    tmp_len_square[(3 * j) + 1] = (point_B[j].width)*(point_B[j].width); 
} 

for (int j = 0; j < N; ++j) { 
    tmp_len_square[(3 * j) + 2] = (point_C[j].height)*(point_C[j].height); 
} 

for (int j = 0; j < N; ++j) { 
    float total = sqrt(tmp_len_square[3 * j] + 
         tmp_len_square[(3 * j) + 1] + 
         tmp_len_square[(3 * j) + 2]); 
    // ... 
} 

Answer 3

您最終可以嘗試優化sqrt函數本身。我建議你看看這個鏈接： Best Square Root Method

Answer 4

將數據重新排列到這一點：

float *pointA_length; 
float *pointB_width; 
float *pointC_height; 

，可能需要你的數據結構的屠宰某種程度，所以你必須選擇不管它是否值得。

現在我們能做的就是這樣寫：

void process_points(float* Alengths, float* Bwidths, float* Cheights, 
        float* output, int n) 
{ 
    for (int i = 0; i < n; i++) { 
     output[i] = sqrt(Alengths[i] * Alengths[i] + 
         Bwidths[i] * Bwidths[i] + 
         Cheights[i] * Cheights[i]); 
    } 
} 

寫像這樣使得它可以自動向量化。例如，針對AVX的GCC和-fno-math-errno -ftree-vectorize可以矢量化該循環。儘管如此，它的確有很多的問題。 __restrict__和對齊屬性只會改善一點。所以這裏有一個手矢量版本，以及：（未測試）

void process_points(float* Alengths, 
        float* Bwidths, 
        float* Cheights, 
        float* output, int n) 
{ 
    for (int i = 0; i < n; i += 8) { 
     __m256 a = _mm256_load_ps(Alengths + i); 
     __m256 b = _mm256_load_ps(Bwidths + i); 
     __m256 c = _mm256_load_ps(Cheights + i); 
     __m256 asq = _mm256_mul_ps(a, a); 
     __m256 sum = _mm256_fmadd_ps(c, c, _mm256_fmadd_ps(b, b, asq)); 
     __m256 hsum = _mm256_mul_ps(sum, _mm256_set1_ps(0.5f)); 
     __m256 invsqrt = _mm256_rsqrt_ps(sum); 
     __m256 s = _mm256_mul_ps(invsqrt, invsqrt); 
     invsqrt = _mm256_mul_ps(sum, _mm256_fnmadd_ps(hsum, s, _mm256_set1_ps(1.5f))); 
     _mm256_store_ps(output + i, _mm256_mul_ps(sum, invsqrt)); 
    } 
} 

這使得一些假設：

• 所有的指針是32對齊。
• n是8的倍數，或者至少緩衝區有足夠的填充，它們永遠不會被超出界限訪問。
• 輸入緩衝區不與輸出緩衝區混淆（它們可能是其中的別名，但是爲什麼）
• 以這種方式計算的平方根的精度稍微降低是可以的（精確到大約22位，而是正確舍入）。
• 與FMADD計算平方的總和可能會稍有不同比如果它使用乘法計算，並補充說，我認爲這沒什麼太
• 目標支持AVX/FMA所以這將實際運行

的方法用於計算這裏使用的平方根是使用近似倒數平方根，改進步驟（y = y * (1.5 - (0.5 * x * y * y))），然後乘以x，因爲x * 1/sqrt(x) = x/sqrt(x) = sqrt(x)。

Answer 5

您的問題可以通過添加更多的上下文來改善。您的代碼是否需要可移植，還是針對特定的編譯器或特定的處理器或處理器系列？也許你願意接受一個通用基線版本，並在運行時選擇特定於目標的優化版本？

此外，您提供的代碼行的上下文很少。它是在一個緊密的循環？還是它散佈在這樣一個循環中的條件代碼中的一堆地方？

我會認爲這是在緊密循環這樣的：

for (int j=0; j<total; ++j) 
    length[j] = sqrt(
     (point_A[j].length)*(point_A[j].length) + 
     (point_B[j].width)*(point_B[j].width) + 
     (point_C[j].height)*(point_C[j].height)); 

我也要去假設你的目標處理器的多核心，該陣列是不同的（或相關元素是不同的），那麼輕鬆取勝是註釋表示OpenMP：

#pragma omp parallel for 
for (int j=0; j<total; ++j) 
    length[j] = sqrt((point_A[j].length)*(point_A[j].length) + 
        (point_B[j].width)*(point_B[j].width) + 
        (point_C[j].height)*(point_C[j].height)); 

編譯g++ -O3 -fopenmp -march=native（或與期望的目標處理器架構替代native）。

如果你知道你的目標，你可能會從gcc標誌-ftree-parallelize-loops=n的並行循環中受益 - 請查看手冊。

現在測量您的績效變化（假設您測量了原始數據，因爲這是一個優化問題）。如果它仍然不夠快，那麼就該考慮更改數據結構，算法或各行代碼。