--- # System prepended metadata title: Parallel Programming and System LAB5 tags: [SP] --- ## CPU 效能調優與演算法設計 > htop在 LAB 4 中顯示了Imbalance的問題，怎麼解? ```c 請輸入 Fibonacci 項數 (建議 30-40): 40 請輸入 Task Cutoff 閾值 (建議 20，輸入 0 代表全平行): 0 [1] 執行純序列計算... 序列結果: 102334155 | 耗時: 0.146667 秒 [2] 執行 OpenMP Task 計算 (Cutoff=0)... 平行結果: 102334155 | 耗時: 215.439894 秒 ``` ## Benchmark : 矩陣乘法 (Naive MatMul) ```c # WSL2 perf 被 ban，所以只能用這個，然後不要跑太大的數字 sudo apt update && sudo apt install valgrind -y valgrind --tool=cachegrind ./matmul_naive # Check export OMP_NUM_THREADS=16 #(default) echo $OMP_NUM_THREADS # EXE gcc fopenmp filename.c -o filename ``` ### Code 2D matrix ```c // matmul_naive_2D.c #include #include #include #include #define N 1000 // 分配不連續的指標陣列 double** alloc_2d(int n) { double **mat = (double **)malloc(n * sizeof(double *)); for (int i = 0; i < n; i++) { mat[i] = (double *)malloc(n * sizeof(double)); } return mat; } void free_2d(double **mat, int n) { for (int i = 0; i < n; i++) free(mat[i]); free(mat); } void init_2d(double **mat, int n) { for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { mat[i][j] = (double)rand() / RAND_MAX; } } } int main() { // 1. 分配分散的記憶體 (Array of Pointers) double **A = alloc_2d(N); double **B = alloc_2d(N); double **C = alloc_2d(N); // C 矩陣需要手動歸零 for(int i=0; i 1D matrix ```c // matmul_naive_1D.c #include #include #include #include #define N 1024 // 矩陣維度 // 初始化矩陣 void init_matrix(double *mat, int size) { for (int i = 0; i < size * size; i++) { mat[i] = (double)rand() / RAND_MAX; } } int main() { // 分配記憶體 (使用一維陣列模擬 N*N) double *A = (double *)malloc(N * N * sizeof(double)); double *B = (double *)malloc(N * N * sizeof(double)); double *C = (double *)calloc(N * N, sizeof(double)); // 初始化為 0 if (!A || !B || !C) { printf("記憶體分配失敗！\n"); return 1; } srand(time(NULL)); init_matrix(A, N); init_matrix(B, N); printf("開始矩陣乘法 (N=%d, Threads=%d)...\n", N, omp_get_max_threads()); double start = omp_get_wtime(); // --- Naive MatMul 核心邏輯 --- #pragma omp parallel for collapse(2) for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { double sum = 0; for (int k = 0; k < N; k++) { // A[i][k] -> A[i*N + k] (橫向，順著走) // B[k][j] -> B[k*N + j] (縱向，跨 N 個元素跳著走) sum += A[i * N + k] * B[k * N + j]; } C[i * N + j] = sum; } } double end = omp_get_wtime(); printf("計算完成。耗時: %.4f 秒\n", end - start); // 簡單校驗 (防止編譯器優化掉整個迴圈) printf("結果範例 C[0]: %f\n", C[0]); free(A); free(B); free(C); return 0; } ``` ```c 1D 版本耗時: 65.4077 秒 ==71886== ==71886== I refs: 26,173,962,133 ==71886== I1 misses: 1,837 ==71886== LLi misses: 1,798 ==71886== I1 miss rate: 0.00% ==71886== LLi miss rate: 0.00% ==71886== ==71886== D refs: 12,074,077,168 (11,055,055,664 rd + 1,019,021,504 wr) ==71886== D1 misses: 1,127,252,106 ( 1,126,001,300 rd + 1,250,806 wr) ==71886== LLd misses: 453,660 ( 77,940 rd + 375,720 wr) ==71886== D1 miss rate: 9.3% ( 10.2% + 0.1% ) ==71886== LLd miss rate: 0.0% ( 0.0% + 0.0% ) ==71886== ==71886== LL refs: 1,127,253,943 ( 1,126,003,137 rd + 1,250,806 wr) ==71886== LL misses: 455,458 ( 79,738 rd + 375,720 wr) ==71886== LL miss rate: 0.0% ( 0.0% + 0.0% ) ``` ### 推測 1. 快取缺失率（Miss Rate）並不正比於最終執行時間 > 2D由於是`malloc`的產物，隨機分布於heap中，相當於前方有做了padding? > 1D則是沒有padding，所以可能導致sharing false? > 如果推論為正確，那padding所下降的miss rate是很可觀的，甚至比是否連續還重要，當然這是對這個scale而言 2. 速度最終還是由1D快了些許，或許表示miss rate所帶來的penalty跟cache conflict比起來還低? 像TLB的cost可能就比想像中的高? 3. `#pragma omp parallel for collapse(2)`的優化 總之，目前的cache 配置導致每讀取一個 double，快取線 (Cache Line) 就會失效一次。所以要透過轉置來解決 ### 轉置 源自於C code 對矩陣記憶體的排列為 Row-major ，使 $$sum += A[i][k] * B[k][j]$$ 這類矩陣乘法在k增加時，B矩陣的存取跨越是一整列，可預測的是，miss rate會降低不少。 >另外，在優化後且不使用compile優化的情況下，16/1 的threads數量並無太大影響 > 額外補充，上述的測試中，LL皆會是0，極有可能是因為L3可以裝入完整的矩陣，故只測量最後一層cache (L3) 的miss rate 為0 ```c // matmul_naive_1TD.c #include #include #include #include #define N 1000 void init_1d(double *mat, int size) { for (int i = 0; i < size * size; i++) { mat[i] = (double)rand() / RAND_MAX; } } int main() { double *A = (double *)malloc(N * N * sizeof(double)); double *B = (double *)malloc(N * N * sizeof(double)); double *B_T = (double *)malloc(N * N * sizeof(double)); // 用於存放轉置後的 B double *C = (double *)calloc(N * N, sizeof(double)); srand(time(NULL)); init_1d(A, N); init_1d(B, N); printf("--- 1D 轉置優化版本 (N=%d) ---\n", N); double start = omp_get_wtime(); // 1. 矩陣轉置階段 (這部分開銷極小，但收益巨大) #pragma omp parallel for collapse(2) for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { B_T[i * N + j] = B[j * N + i]; } } // 2. 核心計算階段 可自行加入num_threads(THREADS)在for後，但測試後16與1的差距並不大 #pragma omp parallel for collapse(2) for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { double sum = 0; for (int k = 0; k < N; k++) { // 現在 A[i*N + k] 和 B_T[j*N + k] 都是連續記憶體存取！ sum += A[i * N + k] * B_T[j * N + k]; } C[i * N + j] = sum; } } double end = omp_get_wtime(); printf("轉置版本耗時: %.4f 秒\n", end - start); free(A); free(B); free(B_T); free(C); return 0; } ``` #### Return ``` 轉置版本耗時: 51.7340 秒 ==75983== ==75983== I refs: 26,197,966,489 ==75983== I1 misses: 1,854 ==75983== LLi misses: 1,815 ==75983== I1 miss rate: 0.00% ==75983== LLi miss rate: 0.00% ==75983== ==75983== D refs: 12,084,079,749 (11,064,056,980 rd + 1,020,022,769 wr) ==75983== D1 misses: 126,629,175 ( 126,128,340 rd + 500,835 wr) ==75983== LLd misses: 627,563 ( 126,821 rd + 500,742 wr) ==75983== D1 miss rate: 1.0% ( 1.1% + 0.0% ) ==75983== LLd miss rate: 0.0% ( 0.0% + 0.0% ) ==75983== ==75983== LL refs: 126,631,029 ( 126,130,194 rd + 500,835 wr) ==75983== LL misses: 629,378 ( 128,636 rd + 500,742 wr) ==75983== LL miss rate: 0.0% ( 0.0% + 0.0% ) ``` ## Tiling > 當矩陣過大，那就切割吧，至少要切到剛剛好L1或L2 !? ### Tiled Matmul Code 即便我們轉置了矩陣，讓存取變得連續，但當 $N$ 變得很大時（例如 $N=2048$ 或更大），整個矩陣的尺寸會遠遠超過你的 **L1/L2 快取**容量。 * **目前的困境**：在計算完矩陣的一列後，當你計算下一列時，之前讀入快取的數據可能已經被「擠掉」了。這意味著同樣的數據會不斷地在 RAM 與 Cache 之間重複搬運。 * **分塊的邏輯**：將大矩陣切成一個個小方塊（Tile），大小剛好能**完全塞進 L1 或 L2 快取**。在處理這個小方塊時，所有的運算都在快取內完成，徹底消除對記憶體頻寬的依賴。 --- 這會將原本的三層迴圈（$i, j, k$）擴展為**六層迴圈**。外三層決定現在處理哪一個「塊」，內三層負責塊內的計算。 > 這裡我們假設分塊大小為 `B` (Block Size)，調整進行測量 ```c #include #include #include #include // 定義矩陣大小與分塊大小 #define N 1000 #define B 8 // Block Size: 根據 L1 Cache 大小調整 (32x32 double 約 8KB) // 輔助函式：取最小值 #define MIN(a, b) ((a) < (b) ? (a) : (b)) void init_matrix(double *mat, int size) { for (int i = 0; i < size * size; i++) { mat[i] = (double)rand() / RAND_MAX; } } int main() { // 1. 記憶體分配 (1D 連續空間) double *A = (double *)malloc(N * N * sizeof(double)); double *B_m = (double *)malloc(N * N * sizeof(double)); // 原始 B double *B_T = (double *)malloc(N * N * sizeof(double)); // 轉置後的 B double *C = (double *)calloc(N * N, sizeof(double)); // 初始化為 0 if (!A || !B_m || !B_T || !C) { printf("記憶體分配失敗！\n"); return -1; } srand(time(NULL)); init_matrix(A, N); init_matrix(B_m, N); printf("--- 最終優化版本: 1D + Transpose + Tiling (N=%d, B=%d) ---\n", N, B); double start_time = omp_get_wtime(); // 2. 矩陣轉置 (將 B 轉置為 B_T，確保計算時為行連續存取) #pragma omp parallel for collapse(2) for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { B_T[i * N + j] = B_m[j * N + i]; } } // 3. 核心計算：六層迴圈 (Tiled Matrix Multiplication) // 外三層迴圈：以「塊」為單位移動 #pragma omp parallel for collapse(2) for (int i = 0; i < N; i += B) { for (int j = 0; j < N; j += B) { for (int k = 0; k < N; k += B) { // 內三層迴圈：在「塊」內部進行計算 for (int ii = i; ii < MIN(i + B, N); ii++) { for (int jj = j; jj < MIN(j + B, N); jj++) { double sum = 0; // 這裡 A 和 B_T 都是 Row-major 連續存取，極利於 Cache L1 與 SIMD for (int kk = k; kk < MIN(k + B, N); kk++) { sum += A[ii * N + kk] * B_T[jj * N + kk]; } C[ii * N + jj] += sum; } } } } } double end_time = omp_get_wtime(); printf("優化版本耗時: %.4f 秒\n", end_time - start_time); // 驗證計算（可選，抽取一個值檢查） // printf("C[500]: %f\n", C[500]); free(A); free(B_m); free(B_T); free(C); return 0; } ``` --- ### 如何決定 Block Size ? 這不是猜出來的，而是根據硬體規格計算出來的。 1. **查詢快取大小**： 在 Linux 下執行 `lscpu` 或 `cat /proc/cpuinfo`。 * 例如：L1 Data Cache = 32 KB。 2. **計算容量限制**： 在計算一個 $B \times B$ 的塊時，我們通常需要同時容納 A 的一塊、B 的一塊（以及 C 的一塊）。 * 如果是 `double` (8 Bytes)，則 $3 \times B^2 \times 8 \text{ Bytes} \leq \text{L1 Cache Size}$。 * 若 L1 為 32 KB，則 $B \approx 32$ 是個不錯的起點。 ### Return ``` // B = 32 with -O3 耗時: 8.7287 秒 ==82241== ==82241== I refs: 5,192,939,004 ==82241== I1 misses: 2,010 ==82241== LLi misses: 1,944 ==82241== I1 miss rate: 0.00% ==82241== LLi miss rate: 0.00% ==82241== ==82241== D refs: 1,132,077,286 (1,082,979,224 rd + 49,098,062 wr) ==82241== D1 misses: 9,969,797 ( 9,593,296 rd + 376,501 wr) ==82241== LLd misses: 628,051 ( 251,857 rd + 376,194 wr) ==82241== D1 miss rate: 0.9% ( 0.9% + 0.8% ) ==82241== LLd miss rate: 0.1% ( 0.0% + 0.8% ) ==82241== ==82241== LL refs: 9,971,807 ( 9,595,306 rd + 376,501 wr) ==82241== LL misses: 629,995 ( 253,801 rd + 376,194 wr) ==82241== LL miss rate: 0.0% ( 0.0% + 0.8% ) ``` ``` // B = 64 without -O3 耗時: 58.0920 秒, (add -O3 耗時: 8.0809 秒, but miss rate increase) ==82745== ==82745== I refs: 31,826,330,518 ==82745== I1 misses: 2,030 ==82745== LLi misses: 1,966 ==82745== I1 miss rate: 0.00% ==82745== LLi miss rate: 0.00% ==82745== ==82745== D refs: 13,348,682,436 (12,283,394,572 rd + 1,065,287,864 wr) ==82745== D1 misses: 7,935,741 ( 7,559,300 rd + 376,441 wr) ==82745== LLd misses: 630,322 ( 254,130 rd + 376,192 wr) ==82745== D1 miss rate: 0.1% ( 0.1% + 0.0% ) ==82745== LLd miss rate: 0.0% ( 0.0% + 0.0% ) ==82745== ==82745== LL refs: 7,937,771 ( 7,561,330 rd + 376,441 wr) ==82745== LL misses: 632,288 ( 256,096 rd + 376,192 wr) ==82745== LL miss rate: 0.0% ( 0.0% + 0.0% ) ``` END LAB 6 https://hackmd.io/@1p0_VOUHTXGcHlBU9A0OEA/rkqFUXtoZx