Samples\0_Introduction\matmul-opt\matmul_opt.cu
/*
* matmul-opt 教学Demo
*
* 本文件包含三种矩阵乘法实现：
* 1. 初学者版本 (Naive) - 直接使用全局内存，简单但性能低
* 2. 优化版本 (Optimized) - 使用共享内存分块+消除Bank冲突+循环展开
* 3. 高级优化版本 (Advanced) - 寄存器分块+float4向量化加载+双缓冲
*
* 目标显卡: NVIDIA Quadro M2000M (sm_50)
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cuda_runtime.h>
#define MATRIX_SIZE 1024
#define CHECK_CUDA_ERROR(err, msg) \
do { \
cudaError_t error = (err); \
if (error != cudaSuccess) { \
printf("CUDA Error at %s: %s (code: %d)\n", msg, cudaGetErrorString(error), error); \
exit(EXIT_FAILURE); \
} \
} while(0)
/*
* CPU版本矩阵乘法：C = A * B
* 这是最基础的三重循环实现，用于生成参考结果
*/
void matrixMulCPU(float* C, const float* A, const float* B,
int M, int N, int K) {
for (int i = 0; i < M; ++i) {
for (int j = 0; j < N; ++j) {
float sum = 0.0f;
for (int k = 0; k < K; ++k) {
sum += A[i * K + k] * B[k * N + j];
}
C[i * N + j] = sum;
}
}
}
/*
* 初学者版本CUDA核函数 - 矩阵乘法
*
* 特点：
* - 每个线程计算C矩阵的一个元素
* - 直接从全局内存读取A和B矩阵
* - 没有使用共享内存
* - 实现简单，但性能低下
*
* 性能瓶颈：
* - 大量重复的全局内存访问
* - 全局内存带宽利用率低
* - 没有数据重用
*/
__global__ void matrixMulNaive(float* C, const float* A, const float* B,
int M, int N, int K) {
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < M && col < N) {
float sum = 0.0f;
for (int k = 0; k < K; ++k) {
sum += A[row * K + k] * B[k * N + col];
}
C[row * N + col] = sum;
}
}
/*
* 优化版本使用的关键参数
* - TILE_DIM: 共享内存分块大小
*/
#define TILE_DIM 32
/*
* 优化版本CUDA核函数 - 矩阵乘法
*
* 优化策略：
* 1. 共享内存分块：将A和B的子矩阵加载到共享内存，减少全局内存访问
* 2. 消除Bank冲突：共享内存声明为 [TILE_DIM][TILE_DIM + 1]
* 3. 循环展开：#pragma unroll 展开内层循环
* 4. 数据复用：每个线程在共享内存中复用数据，避免重复加载
*/
__global__ void matrixMulOptimized(float* C, const float* A, const float* B,
int M, int N, int K) {
__shared__ float As[TILE_DIM][TILE_DIM + 1];
__shared__ float Bs[TILE_DIM][TILE_DIM + 1];
int tx = threadIdx.x;
int ty = threadIdx.y;
int block_row = blockIdx.y * TILE_DIM;
int block_col = blockIdx.x * TILE_DIM;
float sum = 0.0f;
for (int k = 0; k < K; k += TILE_DIM) {
int a_row = block_row + ty;
int a_col = k + tx;
if (a_row < M && a_col < K) {
As[ty][tx] = A[a_row * K + a_col];
} else {
As[ty][tx] = 0.0f;
}
int b_row = k + ty;
int b_col = block_col + tx;
if (b_row < K && b_col < N) {
Bs[ty][tx] = B[b_row * N + b_col];
} else {
Bs[ty][tx] = 0.0f;
}
__syncthreads();
#pragma unroll
for (int e = 0; e < TILE_DIM; ++e) {
sum += As[ty][e] * Bs[e][tx];
}
__syncthreads();
}
int c_row = block_row + ty;
int c_col = block_col + tx;
if (c_row < M && c_col < N) {
C[c_row * N + c_col] = sum;
}
}
/*
* 高级优化版本参数
* 参考: NVIDIA_SGEMM_PRACTICE, LeetCUDA 等开源SGEMM优化项目
*
* - BM, BN: 线程块输出分块大小 (128x128)
* - BK: K维分块大小 (8, 适合float4加载)
* - TM, TN: 每个线程计算的寄存器分块大小 (8x8=64个输出)
* - blockDim = (BM/TM, BN/TN) = (16, 16) = 256 threads
*/
#define BM 128
#define BN 128
#define BK 8
#define TM 8
#define TN 8
/*
* 高级优化版本CUDA核函数 - 矩阵乘法
*
* 三大优化策略：
* 1. 寄存器分块 (Register Tiling): 每个线程计算8x8=64个输出元素
* - 计算访存比从1:2提升到64:16=4:1，大幅减少共享内存访问
* - 中间结果保存在寄存器中，避免反复读写共享内存
*
* 2. float4向量化加载: 每次从全局内存加载4个float (128-bit)
* - 内存事务效率提升4倍
* - 每个线程加载1个float4到A和1个float4到B，恰好覆盖整个分块
*
* 3. 双缓冲 (Double Buffering): 使用2个共享内存缓冲区
* - 当前tile在计算时，下一个tile在加载
* - 重叠计算和内存访问，隐藏全局内存延迟
*/
__global__ void matrixMulAdvanced(float* C, const float* A, const float* B,
int M, int N, int K) {
__shared__ float sA[2][BM][BK + 1];
__shared__ float sB[2][BK][BN + 1];
int tx = threadIdx.x;
int ty = threadIdx.y;
int tid = ty * (BM / TM) + tx;
int bm_start = blockIdx.y * BM;
int bn_start = blockIdx.x * BN;
float rC[TM][TN];
#pragma unroll
for (int i = 0; i < TM; i++) {
#pragma unroll
for (int j = 0; j < TN; j++) {
rC[i][j] = 0.0f;
}
}
int num_k_tiles = (K + BK - 1) / BK;
int a_row_in_tile = tid / 2;
int a_col_in_tile = (tid % 2) * 4;
int b_row_in_tile = tid / 32;
int b_col_in_tile = (tid % 32) * 4;
int a_row_g = bm_start + a_row_in_tile;
int b_col_g = bn_start + b_col_in_tile;
// ---- pre-load tile 0 into buffer 0 ----
int k0 = 0;
int a_col0 = k0 + a_col_in_tile;
if (a_row_g < M && a_col0 < K) {
float4 tmp = *reinterpret_cast<const float4*>(&A[a_row_g * K + a_col0]);
sA[0][a_row_in_tile][a_col_in_tile] = tmp.x;
sA[0][a_row_in_tile][a_col_in_tile + 1] = tmp.y;
sA[0][a_row_in_tile][a_col_in_tile + 2] = tmp.z;
sA[0][a_row_in_tile][a_col_in_tile + 3] = tmp.w;
} else {
sA[0][a_row_in_tile][a_col_in_tile] = 0.0f;
sA[0][a_row_in_tile][a_col_in_tile + 1] = 0.0f;
sA[0][a_row_in_tile][a_col_in_tile + 2] = 0.0f;
sA[0][a_row_in_tile][a_col_in_tile + 3] = 0.0f;
}
int b_row0 = k0 + b_row_in_tile;
if (b_row0 < K && b_col_g < N) {
float4 tmp = *reinterpret_cast<const float4*>(&B[b_row0 * N + b_col_g]);
sB[0][b_row_in_tile][b_col_in_tile] = tmp.x;
sB[0][b_row_in_tile][b_col_in_tile + 1] = tmp.y;
sB[0][b_row_in_tile][b_col_in_tile + 2] = tmp.z;
sB[0][b_row_in_tile][b_col_in_tile + 3] = tmp.w;
} else {
sB[0][b_row_in_tile][b_col_in_tile] = 0.0f;
sB[0][b_row_in_tile][b_col_in_tile + 1] = 0.0f;
sB[0][b_row_in_tile][b_col_in_tile + 2] = 0.0f;
sB[0][b_row_in_tile][b_col_in_tile + 3] = 0.0f;
}
__syncthreads();
// ---- main loop: double buffering ----
for (int kt = 0; kt < num_k_tiles; kt++) {
int buf = kt % 2;
int next_buf = 1 - buf;
if (kt + 1 < num_k_tiles) {
int k_next = (kt + 1) * BK;
int a_col_n = k_next + a_col_in_tile;
if (a_row_g < M && a_col_n < K) {
float4 tmp = *reinterpret_cast<const float4*>(&A[a_row_g * K + a_col_n]);
sA[next_buf][a_row_in_tile][a_col_in_tile] = tmp.x;
sA[next_buf][a_row_in_tile][a_col_in_tile + 1] = tmp.y;
sA[next_buf][a_row_in_tile][a_col_in_tile + 2] = tmp.z;
sA[next_buf][a_row_in_tile][a_col_in_tile + 3] = tmp.w;
} else {
sA[next_buf][a_row_in_tile][a_col_in_tile] = 0.0f;
sA[next_buf][a_row_in_tile][a_col_in_tile + 1] = 0.0f;
sA[next_buf][a_row_in_tile][a_col_in_tile + 2] = 0.0f;
sA[next_buf][a_row_in_tile][a_col_in_tile + 3] = 0.0f;
}
int b_row_n = k_next + b_row_in_tile;
if (b_row_n < K && b_col_g < N) {
float4 tmp = *reinterpret_cast<const float4*>(&B[b_row_n * N + b_col_g]);
sB[next_buf][b_row_in_tile][b_col_in_tile] = tmp.x;
sB[next_buf][b_row_in_tile][b_col_in_tile + 1] = tmp.y;
sB[next_buf][b_row_in_tile][b_col_in_tile + 2] = tmp.z;
sB[next_buf][b_row_in_tile][b_col_in_tile + 3] = tmp.w;
} else {
sB[next_buf][b_row_in_tile][b_col_in_tile] = 0.0f;
sB[next_buf][b_row_in_tile][b_col_in_tile + 1] = 0.0f;
sB[next_buf][b_row_in_tile][b_col_in_tile + 2] = 0.0f;
sB[next_buf][b_row_in_tile][b_col_in_tile + 3] = 0.0f;
}
}
float a_reg[TM];
float b_reg[TN];
#pragma unroll
for (int k = 0; k < BK; k++) {
#pragma unroll
for (int i = 0; i < TM; i++) {
a_reg[i] = sA[buf][ty * TM + i][k];
}
#pragma unroll
for (int j = 0; j < TN; j++) {
b_reg[j] = sB[buf][k][tx * TN + j];
}
#pragma unroll
for (int i = 0; i < TM; i++) {
#pragma unroll
for (int j = 0; j < TN; j++) {
rC[i][j] += a_reg[i] * b_reg[j];
}
}
}
__syncthreads();
}
// ---- write back results ----
#pragma unroll
for (int i = 0; i < TM; i++) {
#pragma unroll
for (int j = 0; j < TN; j++) {
int c_row = bm_start + ty * TM + i;
int c_col = bn_start + tx * TN + j;
if (c_row < M && c_col < N) {
C[c_row * N + c_col] = rC[i][j];
}
}
}
}
/*
* 初始化矩阵
*/
void initMatrix(float* matrix, int size, float value) {
for (int i = 0; i < size; ++i) {
matrix[i] = value;
}
}
/*
* 验证结果正确性
*/
bool verifyResult(const float* gpu_result, const float* cpu_result, int size) {
const float eps = 1e-4f;
int error_count = 0;
for (int i = 0; i < size; ++i) {
float diff = fabs(gpu_result[i] - cpu_result[i]);
if (diff > eps && error_count < 5) {
printf(" Error: idx %d, GPU=%f, CPU=%f, diff=%f\n",
i, gpu_result[i], cpu_result[i], diff);
error_count++;
}
}
if (error_count == 0) {
printf(" Result: PASS\n");
return true;
} else {
printf(" Result: FAIL (showing first %d errors)\n", error_count);
return false;
}
}
/*
* 主函数
*/
int main() {
setvbuf(stdout, NULL, _IONBF, 0);
setvbuf(stderr, NULL, _IONBF, 0);
printf("========================================\n");
printf(" matmul-opt Matrix Multiplication Demo\n");
printf(" GPU: NVIDIA Quadro M2000M (sm_50)\n");
printf("========================================\n\n");
cudaSetDevice(0);
cudaDeviceProp prop;
CHECK_CUDA_ERROR(cudaGetDeviceProperties(&prop, 0), "cudaGetDeviceProperties");
printf("GPU Device: %s\n", prop.name);
printf("Shared Memory Per Block: %zu KB\n", prop.sharedMemPerBlock / 1024);
printf("Warp Size: %d\n", prop.warpSize);
printf("Max Threads Per Block: %d\n", prop.maxThreadsPerBlock);
printf("Multi Processor Count: %d\n\n", prop.multiProcessorCount);
const int M = MATRIX_SIZE;
const int K = MATRIX_SIZE;
const int N = MATRIX_SIZE;
printf("Matrix Size: %dx%d (square)\n", M, K);
printf("Total Elements: %d\n\n", M * N);
float *h_A, *h_B, *h_C_naive, *h_C_opt, *h_C_adv, *h_reference;
CHECK_CUDA_ERROR(cudaMallocHost((void**)&h_A, M * K * sizeof(float)), "cudaMallocHost h_A");
CHECK_CUDA_ERROR(cudaMallocHost((void**)&h_B, K * N * sizeof(float)), "cudaMallocHost h_B");
CHECK_CUDA_ERROR(cudaMallocHost((void**)&h_C_naive, M * N * sizeof(float)), "cudaMallocHost h_C_naive");
CHECK_CUDA_ERROR(cudaMallocHost((void**)&h_C_opt, M * N * sizeof(float)), "cudaMallocHost h_C_opt");
CHECK_CUDA_ERROR(cudaMallocHost((void**)&h_C_adv, M * N * sizeof(float)), "cudaMallocHost h_C_adv");
CHECK_CUDA_ERROR(cudaMallocHost((void**)&h_reference, M * N * sizeof(float)), "cudaMallocHost h_reference");
float *d_A, *d_B, *d_C_naive, *d_C_opt, *d_C_adv;
CHECK_CUDA_ERROR(cudaMalloc((void**)&d_A, M * K * sizeof(float)), "cudaMalloc d_A");
CHECK_CUDA_ERROR(cudaMalloc((void**)&d_B, K * N * sizeof(float)), "cudaMalloc d_B");
CHECK_CUDA_ERROR(cudaMalloc((void**)&d_C_naive, M * N * sizeof(float)), "cudaMalloc d_C_naive");
CHECK_CUDA_ERROR(cudaMalloc((void**)&d_C_opt, M * N * sizeof(float)), "cudaMalloc d_C_opt");
CHECK_CUDA_ERROR(cudaMalloc((void**)&d_C_adv, M * N * sizeof(float)), "cudaMalloc d_C_adv");
printf("Initializing matrices...\n");
initMatrix(h_A, M * K, 1.0f);
initMatrix(h_B, K * N, 0.01f);
printf("Computing CPU reference result...\n");
matrixMulCPU(h_reference, h_A, h_B, M, N, K);
CHECK_CUDA_ERROR(cudaMemcpy(d_A, h_A, M * K * sizeof(float), cudaMemcpyHostToDevice), "cudaMemcpy h_A->d_A");
CHECK_CUDA_ERROR(cudaMemcpy(d_B, h_B, K * N * sizeof(float), cudaMemcpyHostToDevice), "cudaMemcpy h_B->d_B");
dim3 threads_naive(32, 32);
dim3 grid_naive((N + 31) / 32, (M + 31) / 32);
dim3 threads_opt(TILE_DIM, TILE_DIM);
dim3 grid_opt((N + TILE_DIM - 1) / TILE_DIM, (M + TILE_DIM - 1) / TILE_DIM);
dim3 threads_adv(BM / TM, BN / TN);
dim3 grid_adv((N + BN - 1) / BN, (M + BM - 1) / BM);
printf("\nThread Configuration:\n");
printf(" Naive: Block=%dx%d, Grid=%dx%d\n",
threads_naive.x, threads_naive.y, grid_naive.x, grid_naive.y);
printf(" Optim: Block=%dx%d, Grid=%dx%d\n",
threads_opt.x, threads_opt.y, grid_opt.x, grid_opt.y);
printf(" Advanced: Block=%dx%d, Grid=%dx%d\n",
threads_adv.x, threads_adv.y, grid_adv.x, grid_adv.y);
printf("\n");
cudaEvent_t start, stop;
CHECK_CUDA_ERROR(cudaEventCreate(&start), "cudaEventCreate start");
CHECK_CUDA_ERROR(cudaEventCreate(&stop), "cudaEventCreate stop");
const int iterations = 100;
float time_naive = 0.0f, time_opt = 0.0f, time_adv = 0.0f;
printf("========================================\n");
printf(" Running Performance Tests\n");
printf("========================================\n\n");
printf("[Naive Implementation]\n");
matrixMulNaive<<<grid_naive, threads_naive>>>(d_C_naive, d_A, d_B, M, N, K);
CHECK_CUDA_ERROR(cudaGetLastError(), "matrixMulNaive kernel launch");
CHECK_CUDA_ERROR(cudaDeviceSynchronize(), "matrixMulNaive device sync");
printf(" Kernel executed successfully\n");
CHECK_CUDA_ERROR(cudaEventRecord(start), "cudaEventRecord start naive");
for (int i = 0; i < iterations; ++i) {
matrixMulNaive<<<grid_naive, threads_naive>>>(d_C_naive, d_A, d_B, M, N, K);
}
CHECK_CUDA_ERROR(cudaEventRecord(stop), "cudaEventRecord stop naive");
CHECK_CUDA_ERROR(cudaEventSynchronize(stop), "cudaEventSynchronize naive");
CHECK_CUDA_ERROR(cudaEventElapsedTime(&time_naive, start, stop), "cudaEventElapsedTime naive");
CHECK_CUDA_ERROR(cudaMemcpy(h_C_naive, d_C_naive, M * N * sizeof(float), cudaMemcpyDeviceToHost), "cudaMemcpy d_C_naive->h_C_naive");
printf(" First 5 values: %.4f %.4f %.4f %.4f %.4f\n",
h_C_naive[0], h_C_naive[1], h_C_naive[2], h_C_naive[3], h_C_naive[4]);
verifyResult(h_C_naive, h_reference, M * N);
float avg_naive = time_naive / iterations;
double gflops_naive = (2.0 * M * N * K / 1e9) / (avg_naive / 1000.0);
printf(" Average Time: %.4f ms\n", avg_naive);
printf(" Performance: %.2f GFLOPS\n\n", gflops_naive);
printf("[Optimized Implementation]\n");
matrixMulOptimized<<<grid_opt, threads_opt>>>(d_C_opt, d_A, d_B, M, N, K);
CHECK_CUDA_ERROR(cudaGetLastError(), "matrixMulOptimized kernel launch");
CHECK_CUDA_ERROR(cudaDeviceSynchronize(), "matrixMulOptimized device sync");
printf(" Kernel executed successfully\n");
CHECK_CUDA_ERROR(cudaEventRecord(start), "cudaEventRecord start opt");
for (int i = 0; i < iterations; ++i) {
matrixMulOptimized<<<grid_opt, threads_opt>>>(d_C_opt, d_A, d_B, M, N, K);
}
CHECK_CUDA_ERROR(cudaEventRecord(stop), "cudaEventRecord stop opt");
CHECK_CUDA_ERROR(cudaEventSynchronize(stop), "cudaEventSynchronize opt");
CHECK_CUDA_ERROR(cudaEventElapsedTime(&time_opt, start, stop), "cudaEventElapsedTime opt");
CHECK_CUDA_ERROR(cudaMemcpy(h_C_opt, d_C_opt, M * N * sizeof(float), cudaMemcpyDeviceToHost), "cudaMemcpy d_C_opt->h_C_opt");
printf(" First 5 values: %.4f %.4f %.4f %.4f %.4f\n",
h_C_opt[0], h_C_opt[1], h_C_opt[2], h_C_opt[3], h_C_opt[4]);
verifyResult(h_C_opt, h_reference, M * N);
float avg_opt = time_opt / iterations;
double gflops_opt = (2.0 * M * N * K / 1e9) / (avg_opt / 1000.0);
printf(" Average Time: %.4f ms\n", avg_opt);
printf(" Performance: %.2f GFLOPS\n\n", gflops_opt);
printf("[Advanced Implementation]\n");
matrixMulAdvanced<<<grid_adv, threads_adv>>>(d_C_adv, d_A, d_B, M, N, K);
CHECK_CUDA_ERROR(cudaGetLastError(), "matrixMulAdvanced kernel launch");
CHECK_CUDA_ERROR(cudaDeviceSynchronize(), "matrixMulAdvanced device sync");
printf(" Kernel executed successfully\n");
CHECK_CUDA_ERROR(cudaEventRecord(start), "cudaEventRecord start adv");
for (int i = 0; i < iterations; ++i) {
matrixMulAdvanced<<<grid_adv, threads_adv>>>(d_C_adv, d_A, d_B, M, N, K);
}
CHECK_CUDA_ERROR(cudaEventRecord(stop), "cudaEventRecord stop adv");
CHECK_CUDA_ERROR(cudaEventSynchronize(stop), "cudaEventSynchronize adv");
CHECK_CUDA_ERROR(cudaEventElapsedTime(&time_adv, start, stop), "cudaEventElapsedTime adv");
CHECK_CUDA_ERROR(cudaMemcpy(h_C_adv, d_C_adv, M * N * sizeof(float), cudaMemcpyDeviceToHost), "cudaMemcpy d_C_adv->h_C_adv");
printf(" First 5 values: %.4f %.4f %.4f %.4f %.4f\n",
h_C_adv[0], h_C_adv[1], h_C_adv[2], h_C_adv[3], h_C_adv[4]);
verifyResult(h_C_adv, h_reference, M * N);
float avg_adv = time_adv / iterations;
double gflops_adv = (2.0 * M * N * K / 1e9) / (avg_adv / 1000.0);
printf(" Average Time: %.4f ms\n", avg_adv);
printf(" Performance: %.2f GFLOPS\n\n", gflops_adv);
float speedup = time_naive / time_opt;
float speedup_adv = time_naive / time_adv;
printf("========================================\n");
printf(" Performance Summary\n");
printf("========================================\n");
printf(" Naive Time: %.4f ms (%.2f GFLOPS)\n", avg_naive, gflops_naive);
printf(" Optimized Time: %.4f ms (%.2f GFLOPS)\n", avg_opt, gflops_opt);
printf(" Advanced Time: %.4f ms (%.2f GFLOPS)\n", avg_adv, gflops_adv);
printf(" ------------------------------------\n");
printf(" Naive -> Optimized: %.2fx speedup\n", speedup);
printf(" Naive -> Advanced: %.2fx speedup\n", speedup_adv);
printf(" Optimized -> Advanced:%.2fx speedup\n", time_opt / time_adv);
printf("========================================\n\n");
CHECK_CUDA_ERROR(cudaEventDestroy(start), "cudaEventDestroy start");
CHECK_CUDA_ERROR(cudaEventDestroy(stop), "cudaEventDestroy stop");
CHECK_CUDA_ERROR(cudaFreeHost(h_A), "cudaFreeHost h_A");
CHECK_CUDA_ERROR(cudaFreeHost(h_B), "cudaFreeHost h_B");
CHECK_CUDA_ERROR(cudaFreeHost(h_C_naive), "cudaFreeHost h_C_naive");
CHECK_CUDA_ERROR(cudaFreeHost(h_C_opt), "cudaFreeHost h_C_opt");
CHECK_CUDA_ERROR(cudaFreeHost(h_C_adv), "cudaFreeHost h_C_adv");
CHECK_CUDA_ERROR(cudaFreeHost(h_reference), "cudaFreeHost h_reference");
CHECK_CUDA_ERROR(cudaFree(d_A), "cudaFree d_A");
CHECK_CUDA_ERROR(cudaFree(d_B), "cudaFree d_B");
CHECK_CUDA_ERROR(cudaFree(d_C_naive), "cudaFree d_C_naive");
CHECK_CUDA_ERROR(cudaFree(d_C_opt), "cudaFree d_C_opt");
CHECK_CUDA_ERROR(cudaFree(d_C_adv), "cudaFree d_C_adv");
return 0;
}
/*
* ============================================
* 性能优化原理详解
* ============================================
*
* 【初学者版本的性能瓶颈】
*
* 1. 全局内存重复访问：
* - 计算C[i,j]时，需要读取A的第i行和B的第j列
* - 对于K维的每个元素，都要访问全局内存
* - 假设K=256，每个元素计算需要256次全局内存读取
* - 总计：M*N*K次全局内存访问
*
* 2. 内存带宽利用率低：
* - 每个线程独立访问不同的内存地址
* - 同一warp内的线程访问的地址不连续
* - 无法合并内存访问，导致带宽浪费
*
* 【优化版本的核心改进】
*
* 1. 共享内存分块（Tiling）：
* - 将A和B矩阵分割成32x32的小块
* - 一个线程块只处理一个小块
* - 数据加载到共享内存后，在共享内存中完成计算
* - 效果：全局内存访问次数减少为原来的 1/TILE_DIM
* - 原理：共享内存比全局内存快10-100倍
*
* 2. 消除Bank冲突：
* - 共享内存声明 [32][33] 而非 [32][32]
* - 多出的一个float元素作为padding
* - 同一warp的32个线程访问不同的bank
* - 效果：避免共享内存访问冲突，提高带宽
* - 原理：GPU共享内存分为32个bank，同时访问同一bank会串行化
*
* 3. 循环展开：
* - 使用 #pragma unroll 展开内层循环
* - 编译器生成并行的加载和计算指令
* - 效果：减少循环控制开销，提高指令级并行
* - 原理：现代CPU/GPU支持乱序执行和流水线，展开循环可提高效率
*
* 4. 线程块大小优化：
* - 使用 dim3(32, 32) 的线程块，每个线程块刚好一个warp宽
* - 共享内存声明为 [32][33]，消除bank冲突
* - 效果：warp内线程访问共享内存无冲突，提高带宽利用率
* - 原理：32个线程对应32个bank，padding避免列方向冲突
*
* 【为什么这些优化有效？】
*
* 1. 数据重用：
* - 在Naive版本中，A[i,k]被N个线程读取
* - 在优化版本中，A[i,k]只被读取一次到共享内存
* - 然后被TILE_DIM个计算重用
* - 这就是空间局部性的利用
*
* 2. 内存访问模式：
* - Naive版本的内存访问是随机的
* - 优化版本的内存访问是顺序的、合并的
* - 连续的内存访问可以利用缓存和预取
*
* 3. 计算与内存重叠：
* - 加载数据到共享内存和计算可以重叠
* - 一个块在计算时，另一个块可以加载数据
* - 这就是指令级并行的利用
*
* 【理论分析】
*
* 在sm_50架构 (Quadro M2000M) 上：
* - 5个SM (流式多处理器)
* - 每个SM最多支持2048个线程
* - 共享内存带宽：约 15 TB/s
* - 全局内存带宽：约 80-100 GB/s
* - 共享内存比全局内存快 150-187 倍
*
* 通过使用共享内存：
* - 将大部分内存访问从全局内存转移到共享内存
* - 理论上可以获得 10-100 倍的性能提升
*
* 通过消除bank冲突：
* - 共享内存带宽利用率从 1/32 提升到 100%
* - 理论上可以获得 2-32 倍的性能提升
*
* 综合这些优化：
* - 实际性能提升通常在 2-5 倍之间
* - 对于小矩阵(256x256)，提升约1.5-2倍
* - 对于大矩阵(1024x1024)，提升约2-4倍
* - 具体数值取决于矩阵大小和GPU型号
*/
/*
* ============================================
* 高级版本优化原理详解 (matrixMulAdvanced)
* ============================================
*
* 实测性能对比 (1024x1024 矩阵, Quadro M2000M):
* Naive: 46.15 ms ( 46.53 GFLOPS)
* Optimized: 14.33 ms (149.82 GFLOPS)
* Advanced: 3.66 ms (586.72 GFLOPS)
*
* Naive -> Advanced: 12.61x 加速
* Optimized -> Advanced: 3.92x 加速
*
* 【高级版本三大核心优化策略】
*
* 1. 寄存器分块 (Register Tiling)
* - 参数: BM=128, BN=128, BK=8, TM=8, TN=8
* - 每个线程计算 8x8=64 个输出元素 (而非优化版的1个)
* - 线程块大小: (BM/TM) x (BN/TN) = 16x16 = 256 线程
* - 每个线程块输出 128x128=16384 个元素
*
* 为什么能提升性能:
* a) 计算访存比大幅提升:
* - 优化版: 每个输出元素需读2次共享内存 (A和B各1次), 计算访存比 1:2
* - 高级版: 每64个输出元素只需读 8+8=16 次共享内存, 计算访存比 64:16 = 4:1
* - 共享内存访问量减少为优化版的 1/8
* b) 寄存器使用充分:
* - 64个累加器 rC[8][8] 保存在寄存器中 (256 bytes)
* - a_reg[8] 和 b_reg[8] 也在寄存器中 (64 bytes)
* - 中间结果不写回共享内存, 避免反复读写
* c) 指令级并行 (ILP) 更高:
* - 内层循环展开后, 8x8=64 个FMA指令可并行发射
* - 充分利用GPU的指令调度窗口
*
* 2. float4 向量化加载
* - 每次内存事务加载 4 个 float (128-bit = 16 bytes)
* - A tile: 128x8 = 1024 元素 = 256 个 float4, 每个线程加载1个
* - B tile: 8x128 = 1024 元素 = 256 个 float4, 每个线程加载1个
*
* 为什么能提升性能:
* a) 内存事务效率提升4倍:
* - 优化版: 每次加载4字节 (1个float)
* - 高级版: 每次加载16字节 (4个float, 对齐的128-bit事务)
* - GPU全局内存最小事务粒度为32B/128B, float4充分利用带宽
* b) 加载指令数量减少4倍:
* - 相同数据量所需的 LDG 指令减少为1/4
* - 减少指令发射开销和指令缓存压力
* c) 对齐访问:
* - float4 要求16字节对齐, 天然满足合并访问要求
* - 同一warp的32个线程访问连续的128-byte对齐块
*
* 3. 双缓冲 (Double Buffering)
* - 使用两个共享内存缓冲区: sA[2][...], sB[2][...]
* - 前半个循环加载 tile_{t+1} 到 buffer B
* - 后半个循环从 buffer A 计算 tile_t
*
* 为什么能提升性能:
* a) 隐藏全局内存延迟:
* - 优化版: 加载tile -> 等待完成 -> 计算 -> 再加载下一个tile
* - 高级版: 计算当前tile的同时, 后台加载下一个tile
* - 全局内存延迟 (200-400 cycles) 被计算时间掩盖
* b) 流水线化执行:
* - LDG指令 (全局内存加载) 和 FFMA指令 (浮点乘加) 使用不同执行单元
* - 可以真正并行执行, 提高SM利用率
* c) __syncthreads 的巧妙使用:
* - 同步放在计算之后, 确保计算完成后再切换buffer
* - 加载下一个tile不需要额外同步 (在下一次循环迭代开始时自然完成)
*
* 【其他优化细节】
*
* 4. Bank冲突消除 (Padding)
* - sA[BM][BK+1]: 第2维从8改为9, 消除行方向的bank冲突
* - sB[BK][BN+1]: 第2维从128改为129, 消除列方向的bank冲突
* - 8x8分块计算时, 同一warp的32个线程访问8行不同数据, padding确保无冲突
*
* 5. 全循环展开
* - BK=8, TM=8, TN=8 均为编译时常量, #pragma unroll 可完全展开
* - 内层三重循环 8*8*8=512 条指令展开为直线代码
* - 编译器可充分调度指令, 消除分支开销
*
* 6. 分块参数选择依据
* - BM=BN=128: 较大的线程块输出, 增加数据重用, 减少全局内存访问
* - BK=8: 配合float4加载 (8/4=2个float4/行), 每线程正好1个float4
* - TM=TN=8: 每线程64个输出, 平衡寄存器压力与计算访存比
* - 256线程/块: sm_50每SM最多2048线程, 可同时驻留8个块, 充分利用SM
*
* 【理论性能分析】
*
* Quadro M2000M (sm_50) 硬件参数:
* - 5 个 SM, 每SM 128 个 CUDA核心, 共 640 核心
* - 基础频率 ~1072 MHz, 理论峰值: 640 * 1072M * 2 * 4 ≈ 1372 GFLOPS
* - 实测 586.72 GFLOPS, 达到理论峰值的 ~42.7%
*
* 相比优化版的提升来源:
* - 寄存器分块: 计算访存比从1:2提升到4:1 (8倍改善)
* - float4加载: 内存带宽利用率提升约4倍
* - 双缓冲: 全局内存延迟被计算掩盖 (约1.5-2倍改善)
* - 综合效果: 3.92x 加速 (14.33ms -> 3.66ms)
*
* 参考资料:
* - NVIDIA SGEMM Practice (NVIDIA官方SGEMM优化教程)
* - LeetCUDA (https://github.com/DefTruth/CUDA-Learn-Notes)
* - CUDA Matrix Multiplication Optimization (CUTLASS文档)
* - "Programming Massively Parallel Processors" 第10章 矩阵乘法优化
*/
浙公网安备 33010602011771号