爲什麼乘法前的轉置矩陣導致加速很快

Question

我聽說乘法之前的轉置矩陣會因爲緩存局部性而大大提高運算速度。所以我編寫了一個簡單的C++程序，用行主排序來測試它（編譯需要C++ 11和boost）。

結果令人驚訝：7.43秒比0.94秒。但我不明白爲什麼它加快。事實上，在第二個版本（先轉置）中，乘法代碼通過步幅1模式訪問數據，並且具有比第一個更好的局部性。但是，爲了轉置矩陣B，還必須非順序地訪問數據，並導致很多緩存未命中。分配內存和複製數據的開銷也應該是不可忽略的。那麼爲什麼第二個版本會加速代碼呢？

#include <iostream> 
#include <vector> 
#include <boost/timer/timer.hpp> 
#include <random> 

std::vector<int> random_ints(size_t size) 
{ 
    std::vector<int> result; 
    result.reserve(size); 
    std::random_device rd; 
    std::mt19937 engine(rd()); 
    std::uniform_int_distribution<int> dist(0, 100); 
    for (size_t i = 0; i < size; ++i) 
     result.push_back(dist(engine)); 
    return result; 
} 

// matrix A: m x n; matrix B: n x p; matrix C: m x n; 
std::vector<int> matrix_multiply1(const std::vector<int>& A, const std::vector<int>& B, size_t m, size_t n, size_t p) 
{ 
    boost::timer::auto_cpu_timer t; 
    std::vector<int> C(m * p); 
    for (size_t i = 0; i < m; ++i) 
    { 
     for (size_t j = 0; j < p; ++j) 
     { 
      for (size_t k = 0; k < n; ++k) 
      { 
       C[i * m + j] += A[i * m + k] * B[k * n + j]; 
       // B is accessed non-sequentially 
      } 
     } 
    } 
    return C; 
} 

// matrix A: m x n; matrix B: n x p; matrix C: m x n; 
std::vector<int> matrix_multiply2(const std::vector<int>& A, const std::vector<int>& B, size_t m, size_t n, size_t p) 
{ 
    boost::timer::auto_cpu_timer t; 
    std::vector<int> C(m * p), B_transpose(n * p); 

    // transposing B 
    for (size_t i = 0; i < n; ++i) 
    { 
     for (size_t j = 0; j < p; ++j) 
     { 
      B_transpose[i + j * p] = B[i * n + j]; 
      // B_transpose is accessed non-sequentially 
     } 
    } 

    // multiplication 
    for (size_t i = 0; i < m; ++i) 
    { 
     for (size_t j = 0; j < p; ++j) 
     { 
      for (size_t k = 0; k < n; ++k) 
      { 
       C[i * m + j] += A[i * m + k] * B_transpose[k + j * p]; 
       // all sequential access 
      } 
     } 
    } 
    return C; 
} 

int main() 
{ 
    const size_t size = 1 << 10; 
    auto A = random_ints(size * size); 
    auto C = matrix_multiply1(A, A, size, size, size); 
    std::cout << C.front() << ' ' << C.back() << std::endl; // output part of the result 
    C = matrix_multiply2(A, A, size, size, size); 
    std::cout << C.front() << ' ' << C.back() << std::endl; // compare with output of algorithm 1 
    return 0; 
} 

Answer 1

乘法涉及比轉置更多的訪問，所以它支配執行時間。

您可以只用在尋找環頭，很清楚地看到這一點：

// transpose 
for (size_t i = 0; i < n; ++i) 
    for (size_t j = 0; j < p; ++j) 
     ... 

// multiplication 
for (size_t i = 0; i < m; ++i) 
    for (size_t j = 0; j < p; ++j) 
     for (size_t k = 0; k < n; ++k) 
      ... 

有了一個額外的嵌套，二是要明顯得多的工作。