優化四重嵌套「for」循環

Question

我正在開發一個2D數值模型，並且我想加快特定成員函數的速度，這會減慢我的代碼。該函數需要循環遍歷模型中的每個網格點，然後在每個網格點上執行雙重求和，範圍爲l和m。功能如下：

int Class::Function(void) { 
    double loadingEta; 
    int i,j,l,m; 

    //etaLatLen=64, etaLonLen=2*64 
    //l_max = 12 

    for (i=0; i<etaLatLen; i++) { 
     for (j=0; j < etaLonLen; j++) { 
      loadingEta = 0.0; 
      for (l=0; l<l_max+1; l++) { 
       for (m=0; m<=l; m++) { 
        loadingEta += etaLegendreArray[i][l][m] * (SH_C[l][m]*etaCosMLon[j][m] + SH_S[l][m]*etaSinMLon[j][m]); 
       } 
      } 
      etaNewArray[i][j] = loadingEta; 
     } 
    } 

    return 1; 
} 

我一直在嘗試改變循環順序來加快速度，但無濟於事。任何幫助將非常感激。謝謝！

編輯1：

所有五個數組在我的類的構造函數分配如下：

etaLegendreArray = new double**[etaLatLen]; 
for (int i=0; i<etaLatLen; i++) { 
    etaLegendreArray[i] = new double*[l_max+1]; 
    for (int l=0; l<l_max+1; l++) { 
     etaLegendreArray[i][l] = new double[l_max+1]; 
    } 
} 

SH_C = new double*[l_max+1]; 
SH_S = new double*[l_max+1]; 
for (int i=0; i<l_max+1; i++) { 
    SH_C[i] = new double[l_max+1]; 
    SH_S[i] = new double[l_max+1]; 
} 

etaCosMLon = new double*[etaLonLen]; 
etaSinMLon = new double*[etaLonLen]; 
for (int j=0; j<etaLonLen; j++) { 
    etaCosMLon[j] = new double[l_max+1]; 
    etaSinMLon[j] = new double[l_max+1]; 
} 

也許會更好，如果這些是一維數組，而不是多維？

Answer 1

這裏跳到X-Y領土。我們嘗試加快數據訪問速度，而不是加速算法。

etaLegendreArray = new double**[etaLatLen]; 
for (int i=0; i<etaLatLen; i++) { 
    etaLegendreArray[i] = new double*[l_max+1]; 
    for (int l=0; l<l_max+1; l++) { 
     etaLegendreArray[i][l] = new double[l_max+1]; 
    } 
} 

不創建3D陣列double s。它爲指向數組double的指針數組創建一個指針數組。每個數組都是它自己的內存塊，誰知道它將在哪裏存儲。這導致了一個稱爲「poor spacial locality」的數據結構。所有的結構件可能散落在各處。在三維陣列中，您可以跳到三個不同的地方，以查明您的價值在哪裏。

由於模擬3D陣列所需的許多存儲塊可能遠不及彼此，因此CPU可能無法提前有效加載高速緩存（高速存儲器），必須停止有用的工作它正在等待訪問速度較慢的存儲，可能更頻繁地訪問RAM。這是一個很好的，高水平的表現。另一方面，如果整個數組位於一個內存塊中，是「連續的」，則CPU可以讀取更大的內存塊，也許是所有內存塊，它需要立即進入緩存。此外，如果編譯器知道程序將使用的內存全部在一個大塊中，它可以執行各種常規優化，這將使您的程序更快。

那麼，我們如何獲得一個3D存儲器塊？如果大小是靜態的，這很容易

double etaLegendreArray[SIZE1][SIZE2][SIZE3]; 

這不看是你的情況，所以你想做的事是分配一維數組，因爲這將是一個內存連續塊。

double * etaLegendreArray= new double [SIZE1*SIZE2*SIZE3]; 

，做手工的數組索引數學

etaLegendreArray[(x * SIZE2 + y) * SIZE3 + z] = data; 

貌似這應該是所有額外的數學比較慢，對吧？原來編譯器隱藏的數學看起來很像你每次使用[]時的情況。你幾乎沒有損失，當然不會像失去一個不必要的cache miss那樣多。

但是，在整個地方重複這個數學是瘋狂的，遲早你會搞砸了，即使可讀性耗盡並不是你首先希望死亡的，所以你真的想把1D數組包裝在一個班級幫手處理你的數學。一旦你這樣做了，你可能會讓這個類處理分配和釋放，這樣你就可以利用all that RAII goodness。沒有更多for環路new s和delete s遍佈各處。它全部包裹起來並用弓綁起來。

Here is an example of a 2D Matrix class easily extendable to 3D.這將以一個很好的可預測和緩存友好的方式照顧您可能需要的基本功能。

Answer 2

如果CPU支持它並且編譯器進行了足夠優化，您可能會從the C99 fma（融合乘加）函數中獲得一些小的增益，將一些步驟的操作（乘，然後加）轉換爲一步操作。這也可以提高準確性，因爲對於融合操作你只進行一次浮點舍入，而不是一次乘法和一次加法。

loadingEta += etaLegendreArray[i][l][m] * (SH_C[l][m]*etaCosMLon[j][m] + SH_S[l][m]*etaSinMLon[j][m]); 

到（注意沒有用+=現在，它在fma真實於此）：

loadingEta = fma(etaLegendreArray[i][l][m], fma(SH_C[l][m], etaCosMLon[j][m], SH_S[l][m]*etaSinMLon[j][m]), loadingEta); 

假設我讀它的權利，你可以從改變你的內部循環的表情我不希望任何神奇的性能方面的東西，但它可能會有所幫助（再一次，只有優化足夠讓編譯器內聯硬件指令來完成這項工作;如果它調用一個庫函數，您將失去任何改進到函數調用開銷）。再次，它應該通過避免兩個四捨五入的步驟來提高準確性。

請注意，在some compilers with appropriate compilation flags, they'll convert your original code to hardware FMA instructions for you;如果這是一個選項，我會去，因爲（如你所見）fma功能傾向於減少代碼的可讀性。

您的編譯器也可能提供浮點指令的矢量化版本，這可能會顯着提高性能（請參閱上一個鏈接自動轉換爲FMA）。

大多數其他改進將需要更多關於目標，使用的輸入數組的性質等信息。簡單的線程可能會帶給你一些東西，OpenMP編譯指示可能是一種可以簡化並行化循環的方法（ S）。