一、矩陣乘法數學原理與性能瓶頸

1.1 數學原理

矩陣乘法定義為：給定兩個矩陣 $\mathrm{A}（m\times n）$和$\mathrm{B}（n\times p）$，它們的乘積 $\mathrm{C}=A\times B$ 是一個 $\mathrm{m}\times p$ 的矩陣，其中：

$$C_{i,j}=\sum _{k=1}^n{A}_{i,k}?B_{k,j}$$

每個元素$C[i][j]$ 需要 $n$ 次乘法和 $n?1$ 次加法，總時間復雜度為 $O(mnp)$ 。對于兩個$n\times n$ 方陣，時間復雜度為 $O(n^3)$。

2.2 ?性能問題

2.2.1 內存訪問效率

• ?緩存未命中：大矩陣無法完全放入緩存，導致頻繁訪問內存，增加延遲。
• 非連續訪問：訪問$B$的列時，若存儲為行主序，會導致緩存不友好。

2.2.2 并行化效率

• ?任務分配不均：靜態分配可能導致負載不均衡。
• 同步開銷：多線程間同步可能引入額外開銷。

2.2.3 指令級并行與向量化

• SIMD利用率低：未充分利用向量指令（如AVX、NEON）。

2.2.4 循環結構與數據布局

• 低效循環順序：傳統三重循環（i-j-k）導致內存訪問不連續。

2.3 優化策略

1. 分塊處理：將矩陣劃分為適合緩存的小塊，減少內存訪問。
2. 多線程并行：使用OpenMP、pthread等多線程庫并行計算。
3. SIMD向量化：利用AVX/SSE指令加速計算。
4. 數據布局優化：轉置矩陣 $B$ 或調整存儲順序（行/列主序）。

---

二、分塊策略優化（Cache-Friendly）?

分塊（Tiling）策略是優化矩陣乘法性能的核心技術之一，其核心思想是通過**數據局部性（Locality）**?提升緩存利用率，減少內存訪問延遲。

2.1 為什么需要分塊？

矩陣乘法的樸素實現（三重循環）通常存在以下問題：

• 內存訪問模式差：對矩陣$B$的列訪問（行主序存儲下）不連續，導致**緩存行（Cache Line）**利用率低。
• ?數據重用率低：每個元素$A[i][k]$和$B[k][j]$僅被使用一次，無法利用緩存的時間局部性（Temporal Locality）。
• 緩存容量不足：當矩陣大小超過緩存容量時，頻繁的緩存未命中（Cache Miss）會觸發內存訪問，導致性能驟降。
  分塊策略通過將大矩陣拆分為適合緩存的小塊，使得每個小塊的數據能長時間駐留在緩存中，從而提升數據重用率。

2.2 分塊策略的原理

2.2.1 數據局部性優化

• 時間局部性：同一塊內的數據被多次使用（例如，分塊后$A$的子塊行與$B$的子塊列多次參與計算）。
• ?空間局部性：連續訪問內存中的數據（例如，按行優先順序訪問子塊內的元素）。

2.2.2 分塊步驟

1. 劃分矩陣：將矩陣$A$和$B$劃分為大小為$T\times T$的子塊（Block），$A$按行分塊，$B$按列分塊，$C$對應分塊，分塊大小$T$需根據緩存容量（如L1/L2 Cache大小）調整。
2. 分塊乘法：對每個子塊執行矩陣乘法（子塊乘積累加到$C$的對應位置）。
3. 循環順序調整：將循環順序從樸素的 ?i-j-k 調整為 ?i-k-j，優先處理分塊內的計算。

2.2.3 分塊后的計算流程

for i in 0 to m step T:           # 分塊行循環for j in 0 to p step T:       # 分塊列循環for k in 0 to n step T:   # 分塊中間循環# 計算子塊 C[i:i+T, j:j+T] += A[i:i+T, k:k+T] * B[k:k+T, j:j+T]for ii in i to i+T:   # 子塊內行循環for kk in k to k+T:for jj in j to j+T:C[ii][jj] += A[ii][kk] * B[kk][jj]

2.3 示例代碼

#include <immintrin.h>void block_matmul(int m, int n, int p, float* A, float* B, float* C) {const int T = 64; // 分塊大小for (int i = 0; i < m; i += T) {for (int j = 0; j < p; j += T) {for (int k = 0; k < n; k += T) {// 子塊范圍int i_end = std::min(i + T, m);int j_end = std::min(j + T, p);int k_end = std::min(k + T, n);// 子塊內計算（使用AVX2）for (int ii = i; ii < i_end; ii++) {for (int kk = k; kk < k_end; kk++) {__m256 a = _mm256_broadcast_ss(&A[ii * n + kk]);for (int jj = j; jj < j_end; jj += 8) {__m256 b = _mm256_loadu_ps(&B[kk * p + jj]);__m256 c = _mm256_loadu_ps(&C[ii * p + jj]);c = _mm256_fmadd_ps(a, b, c);_mm256_storeu_ps(&C[ii * p + jj], c);}}}}}}
}

---

三、矩陣轉置

矩陣轉置是線性代數中的基本操作，指將矩陣的行和列互換。具體來說，對于一個$m\times n$的矩陣 $A$，其轉置矩陣 $A^T$是一個$n\times m$的矩陣。

3.1 8x8塊轉置代碼實現

void transpose_8x8_block(size_t i, size_t j, Matrix& result) const {const float* src = data_ + i * cols_ + j;__m256 row0 = _mm256_loadu_ps(src + 0 * cols_);__m256 row1 = _mm256_loadu_ps(src + 1 * cols_);__m256 row2 = _mm256_loadu_ps(src + 2 * cols_);__m256 row3 = _mm256_loadu_ps(src + 3 * cols_);__m256 row4 = _mm256_loadu_ps(src + 4 * cols_);__m256 row5 = _mm256_loadu_ps(src + 5 * cols_);__m256 row6 = _mm256_loadu_ps