爲什麼轉置512x512的矩陣要比轉置513x513的矩陣慢得多？

Question

在不同尺寸的矩陣上進行了一些實驗之後，出現了一種模式。總體上，轉置大小爲2^n的矩陣比轉置大小2^n+1中的一個要慢。對於n的小數值，差異不是很大。然而

大的差異發生在512的值（至少對我來說）

免責聲明：我知道這個功能實際上並沒有轉，因爲元素的雙交換矩陣，但它沒有區別。

下面的代碼：

#define SAMPLES 1000 
#define MATSIZE 512 

#include <time.h> 
#include <iostream> 
int mat[MATSIZE][MATSIZE]; 

void transpose() 
{ 
    for (int i = 0 ; i < MATSIZE ; i++) 
    for (int j = 0 ; j < MATSIZE ; j++) 
    { 
     int aux = mat[i][j]; 
     mat[i][j] = mat[j][i]; 
     mat[j][i] = aux; 
    } 
} 

int main() 
{ 
    //initialize matrix 
    for (int i = 0 ; i < MATSIZE ; i++) 
    for (int j = 0 ; j < MATSIZE ; j++) 
     mat[i][j] = i+j; 

    int t = clock(); 
    for (int i = 0 ; i < SAMPLES ; i++) 
     transpose(); 
    int elapsed = clock() - t; 

    std::cout << "Average for a matrix of " << MATSIZE << ": " << elapsed/SAMPLES; 
} 

更改MATSIZE讓我們改變大小（廢話！）。我張貼在ideone兩個版本：

• 尺寸512 - 平均2.46毫秒 - http://ideone.com/1PV7m
• 尺寸513 - 平均0.75毫秒 - http://ideone.com/NShpo

在我的環境（ MSVS 2010，全面優化），差別類似：

• 尺寸512 - 平均2.19毫秒
• 尺寸513 - 平均0.57毫秒

這究竟是爲什麼？

Answer 1

這個解釋來自於Optimizing software in C++中的Agner Fog，它減少了數據如何被訪問和存儲在緩存中。

有關條款和詳細信息，請參閱wiki entry on caching，我將在此處縮小範圍。

高速緩存組織在集和行。一次只能使用一套，其中包含的任何一行都可以使用。一行可以鏡像的行數乘以行數給我們緩存大小。

對於特定的存儲器地址，就可以計算出該設定它應該與式進行鏡像中：

set = (address/lineSize) % numberOfsets 

這類式是使整個組理想地均勻分佈，因爲每個存儲器地址是作爲很可能會被閱讀（我說理想情況下）。

很明顯，重疊可能發生。在高速緩存未命中的情況下，將在高速緩存中讀取內存並替換舊值。記住每個集合都有許多行，最近最少使用的行將被新讀取的內存覆蓋。

我將設法有所遵循從昂納的例子：

假定每個組有4行，每行保持64個字節。我們首先嚐試讀取地址0x2710，該地址在集合28中。然後我們也嘗試讀取地址0x2F00，0x3700，0x3F00和​​。所有這些屬於同一組。在閱讀​​之前，該集合中的所有行將被佔用。讀取該內存會清除集合中現有的一行，該行最初持有0x2710。問題在於我們讀取的地址是（此例）0x800。這是關鍵步幅（再次，對於這個例子）。

臨界步幅也可以計算：

criticaStride = numberOfSets * lineSize 

變量隔開criticalStride或多個分開爭用相同的高速緩存行。

這是理論部分。接下來，解釋（也是Agner，我正在密切關注以避免犯錯）：

假設一個矩陣爲64x64（記住，效果因緩存而異），一個8kb緩存，每組4行*行大小爲64字節。每行可以容納矩陣中的8個元素（64位int）。

關鍵跨步將是2048字節，這對應於矩陣的4行（在內存中是連續的）。

假設我們正在處理第28行。我們試圖獲取該行的元素，並將它們與第28列的元素交換。行的前8個元素組成緩存行，但它們會進入第28列中的8個不同的緩存行。請記住，關鍵步幅相隔4行（一列中有4個連續元素）。

當在列中到達元素16時（每組4個高速緩存行&間隔4行=故障），ex-0元素將從高速緩存中逐出。當我們到達列的末尾時，所有先前的緩存行將會丟失，並且在訪問下一個元素時需要重新加載（整行被覆蓋）。

其尺寸不是關鍵的步幅的倍數攪亂這個完美方案災難，因爲我們不再使用，除了是至關重要的步幅上垂直件處理，所以重新加載緩存的數量嚴重減少。

另一個免責聲明 - 我只是對解釋有所瞭解，並希望我能指出它，但我可能會誤會。無論如何，我正在等待Mysticial的回覆（或確認）。 :)

Answer 2

Luchian給出解釋爲什麼這種行爲發生，但我認爲這會是一個不錯的主意，以顯示一種可能的解決了這個問題，並在同一時間出示了一下有關緩存忘卻的算法。

你的算法基本上沒有：

for (int i = 0; i < N; i++) 
    for (int j = 0; j < N; j++) 
     A[j][i] = A[i][j]; 

這太可怕了現代CPU。一個解決方案是知道你的緩存系統的細節，並調整算法以避免這些問題。只要你知道那些細節，工作很棒..不是特別便攜。

我們可以做得更好嗎？是的，我們可以：這個問題的一般方法是cache oblivious algorithms，作爲它的名字說避免依賴於特定的緩存大小[1]

該解決方案是這樣的：

void recursiveTranspose(int i0, int i1, int j0, int j1) { 
    int di = i1 - i0, dj = j1 - j0; 
    const int LEAFSIZE = 32; // well ok caching still affects this one here 
    if (di >= dj && di > LEAFSIZE) { 
     int im = (i0 + i1)/2; 
     recursiveTranspose(i0, im, j0, j1); 
     recursiveTranspose(im, i1, j0, j1); 
    } else if (dj > LEAFSIZE) { 
     int jm = (j0 + j1)/2; 
     recursiveTranspose(i0, i1, j0, jm); 
     recursiveTranspose(i0, i1, jm, j1); 
    } else { 
    for (int i = i0; i < i1; i++) 
     for (int j = j0; j < j1; j++) 
      mat[j][i] = mat[i][j]; 
    } 
} 

稍微複雜一些，但一個簡短的測試顯示MATSIZE 8192

int main() { 
    LARGE_INTEGER start, end, freq; 
    QueryPerformanceFrequency(&freq); 
    QueryPerformanceCounter(&start); 
    recursiveTranspose(0, MATSIZE, 0, MATSIZE); 
    QueryPerformanceCounter(&end); 
    printf("recursive: %.2fms\n", (end.QuadPart - start.QuadPart)/(double(freq.QuadPart)/1000)); 

    QueryPerformanceCounter(&start); 
    transpose(); 
    QueryPerformanceCounter(&end); 
    printf("iterative: %.2fms\n", (end.QuadPart - start.QuadPart)/(double(freq.QuadPart)/1000)); 
    return 0; 
} 

results: 
recursive: 480.58ms 
iterative: 3678.46ms 

的東西在我的古E8400與VS2010發佈的x64挺有意思的，testcode編輯：關於大小的影響：它是那麼明顯，雖然仍明顯在一定程度上，這是因爲我們將迭代解決方案用作葉節點，而不是遞歸到1（通常遞歸算法的優化）。如果我們設置LEAFSIZE = 1，緩存對我沒有影響[8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms - 這是在誤差範圍內，波動在100ms區域;如果我們想要完全準確的數值，這個「基準」並不是我會感到太舒服的原因]）

[1]這個東西的來源：好吧，如果你不能從一個合作過的人Leiserson和co ..我認爲他們的論文是一個很好的起點。這些算法仍然很少被描述 - CLR有一個關於它們的腳註。儘管如此，這仍然是給人們驚喜的好方法。

---

編輯（注：我不是誰張貼了這個答案的一個;我只是想補充這一點）：
這裏是上面代碼的完整C++版本：

template<class InIt, class OutIt> 
void transpose(InIt const input, OutIt const output, 
    size_t const rows, size_t const columns, 
    size_t const r1 = 0, size_t const c1 = 0, 
    size_t r2 = ~(size_t) 0, size_t c2 = ~(size_t) 0, 
    size_t const leaf = 0x20) 
{ 
    if (!~c2) { c2 = columns - c1; } 
    if (!~r2) { r2 = rows - r1; } 
    size_t const di = r2 - r1, dj = c2 - c1; 
    if (di >= dj && di > leaf) 
    { 
     transpose(input, output, rows, columns, r1, c1, (r1 + r2)/2, c2); 
     transpose(input, output, rows, columns, (r1 + r2)/2, c1, r2, c2); 
    } 
    else if (dj > leaf) 
    { 
     transpose(input, output, rows, columns, r1, c1, r2, (c1 + c2)/2); 
     transpose(input, output, rows, columns, r1, (c1 + c2)/2, r2, c2); 
    } 
    else 
    { 
     for (ptrdiff_t i1 = (ptrdiff_t) r1, i2 = (ptrdiff_t) (i1 * columns); 
      i1 < (ptrdiff_t) r2; ++i1, i2 += (ptrdiff_t) columns) 
     { 
      for (ptrdiff_t j1 = (ptrdiff_t) c1, j2 = (ptrdiff_t) (j1 * rows); 
       j1 < (ptrdiff_t) c2; ++j1, j2 += (ptrdiff_t) rows) 
      { 
       output[j2 + i1] = input[i2 + j1]; 
      } 
     } 
    } 
} 

Answer 3

作爲對Luchian Grigore's answer中解釋的說明，以下是64x64和65x65矩陣這兩種情況下的矩陣緩存存在情況（請參閱上面的鏈接，瞭解有關數字的詳細信息）。下面

色彩的動畫含義如下：

• - 不在緩存中，
• - 在高速緩存中，
• - 高速緩存命中，
• - 從RAM剛讀，
• - 緩存未命中。

64×64的情況下：

注意如何幾乎每一個訪問緩存未命中的一個新行的結果。現在怎麼它看起來正常情況下，一個65x65矩陣：

在這裏你可以看到，大部分的初始磨合後訪問的高速緩存命中。這就是CPU緩存一般用於如何工作的方式。