当前位置：首页 > news >正文

(done) CUDA 和 CPU 性能对比，矩阵加法和矩阵乘法对比

news 2025/9/12 15:12:47

矩阵加法 cuda/CPU 对比

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <cuda_runtime.h>// CUDA 矩阵加法内核函数
__global__ void matrixAddCUDA(int *a, int *b, int *c, int width, int height)
{int col = blockIdx.x * blockDim.x + threadIdx.x;int row = blockIdx.y * blockDim.y + threadIdx.y;if (col < width && row < height) {int index = row * width + col;c[index] = a[index] + b[index];}
}// CPU 矩阵加法函数
void matrixAddCPU(int *a, int *b, int *c, int width, int height)
{for (int row = 0; row < height; row++) {for (int col = 0; col < width; col++) {int index = row * width + col;c[index] = a[index] + b[index];}}
}// 初始化矩阵
void initializeMatrix(int *matrix, int size)
{for (int i = 0; i < size; i++) {matrix[i] = rand() % 100; // 0-99 的随机数}
}// 验证结果是否正确
bool verifyResult(int *cpu_result, int *gpu_result, int size)
{for (int i = 0; i < size; i++) {if (cpu_result[i] != gpu_result[i]) {printf("验证失败! 索引 %d: CPU=%d, GPU=%d\n", i, cpu_result[i], gpu_result[i]);return false;}}return true;
}int main()
{// 设置矩阵大小const int WIDTH = 1024;const int HEIGHT = 1024;const int SIZE = WIDTH * HEIGHT;const int BYTES = SIZE * sizeof(int); // 字节数量printf("矩阵大小: %d x %d = %d 个元素\n", WIDTH, HEIGHT, SIZE);printf("总数据量: %.2f MB\n", BYTES / (1024.0 * 1024.0) * 2); // 两个输入矩阵// 分配主机内存int *h_a = (int*)malloc(BYTES); // 输入矩阵 Aint *h_b = (int*)malloc(BYTES); // 输入矩阵 Bint *h_c_cpu = (int*)malloc(BYTES); // CPU 结果矩阵int *h_c_gpu = (int*)malloc(BYTES); // GPU 结果矩阵// 随机初始化输入矩阵srand(time(NULL));initializeMatrix(h_a, SIZE);initializeMatrix(h_b, SIZE);// 分配GPU设备内存int *d_a, *d_b, *d_c;cudaMalloc((void**)&d_a, BYTES);cudaMalloc((void**)&d_b, BYTES);cudaMalloc((void**)&d_c, BYTES);// 拷贝数据到设备cudaMemcpy(d_a, h_a, BYTES, cudaMemcpyHostToDevice);cudaMemcpy(d_b, h_b, BYTES, cudaMemcpyHostToDevice);// 设置 CUDA 线程块和网格维度dim3 blockSize(16, 16); // 256 个线程每块dim3 gridSize((WIDTH + blockSize.x - 1) / blockSize.x, (HEIGHT + blockSize.y - 1) / blockSize.y);printf("\nCUDA 配置: 网格(%d, %d), 块(%d, %d)\n", gridSize.x, gridSize.y, blockSize.x, blockSize.y);// 测试 CPU 版本 (生成参考结果)clock_t cpu_start = clock();matrixAddCPU(h_a, h_b, h_c_cpu, WIDTH, HEIGHT);clock_t cpu_end = clock();double cpu_time = ((double)(cpu_end - cpu_start)) / CLOCKS_PER_SEC * 1000.0;// 测试 GPU 版本cudaEvent_t start, stop;cudaEventCreate(&start);cudaEventCreate(&stop);cudaEventRecord(start);// 启动 CUDA 内核 (运行 kernel)matrixAddCUDA<<<gridSize, blockSize>>>(d_a, d_b, d_c, WIDTH, HEIGHT);cudaEventRecord(stop);cudaEventSynchronize(stop);// 拷贝结果回主机cudaMemcpy(h_c_gpu, d_c, BYTES, cudaMemcpyDeviceToHost);// 计算 GPU 执行时间float gpu_time = 0;cudaEventElapsedTime(&gpu_time, start, stop);// 验证结果bool isCorrect = verifyResult(h_c_cpu, h_c_gpu, SIZE);// 输出性能结果printf("\n性能比较:\n");printf("CPU 执行时间: %.2f ms\n", cpu_time);printf("GPU 执行时间: %.2f ms\n", gpu_time);printf("加速比: %.2f 倍\n", cpu_time / gpu_time);printf("结果验证: %s\n", isCorrect ? "正确" : "错误");// 清理资源free(h_a);free(h_b);free(h_c_cpu);free(h_c_gpu);cudaFree(d_a);cudaFree(d_b);cudaFree(d_c);cudaEventDestroy(start);cudaEventDestroy(stop);return 0;
}

矩阵乘法 cuda/CPU 对比

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <cuda_runtime.h>// 定义矩阵大小（可以调整以测试不同规模）
#define N 1024
#define BLOCK_SIZE 16// CPU矩阵乘法实现
void matrixMulCPU(float *A, float *B, float *C, int n) {for (int i = 0; i < n; i++) {for (int j = 0; j < n; j++) {float sum = 0.0f;for (int k = 0; k < n; k++) {sum += A[i * n + k] * B[k * n + j];}C[i * n + j] = sum;}}
}// 简单的CUDA矩阵乘法核函数（不使用共享内存）
__global__ void matrixMulSimple(float *A, float *B, float *C, int n) {int row = blockIdx.y * blockDim.y + threadIdx.y;int col = blockIdx.x * blockDim.x + threadIdx.x;if (row < n && col < n) {float sum = 0.0f;for (int k = 0; k < n; k++) {sum += A[row * n + k] * B[k * n + col];}C[row * n + col] = sum;}
}// 使用共享内存优化的CUDA矩阵乘法核函数
__global__ void matrixMulShared(float *A, float *B, float *C, int n) {// 为当前块分配共享内存__shared__ float As[BLOCK_SIZE][BLOCK_SIZE];__shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];// 计算线程的行列索引int row = blockIdx.y * BLOCK_SIZE + threadIdx.y;int col = blockIdx.x * BLOCK_SIZE + threadIdx.x;float sum = 0.0f;// 循环遍历所有子矩阵for (int k = 0; k < n; k += BLOCK_SIZE) {// 将数据加载到共享内存if (row < n && (k + threadIdx.x) < n) {As[threadIdx.y][threadIdx.x] = A[row * n + k + threadIdx.x];} else {As[threadIdx.y][threadIdx.x] = 0.0f;}if ((k + threadIdx.y) < n && col < n) {Bs[threadIdx.y][threadIdx.x] = B[(k + threadIdx.y) * n + col];} else {Bs[threadIdx.y][threadIdx.x] = 0.0f;}// 等待所有线程完成数据加载__syncthreads();// 计算当前子矩阵的乘积for (int i = 0; i < BLOCK_SIZE; i++) {sum += As[threadIdx.y][i] * Bs[i][threadIdx.x];}// 等待所有线程完成计算__syncthreads();}// 将结果写入全局内存if (row < n && col < n) {C[row * n + col] = sum;}
}// 初始化矩阵
void initMatrix(float *matrix, int n, int seed) {srand(seed);for (int i = 0; i < n * n; i++) {matrix[i] = (float)rand() / RAND_MAX;}
}// 验证结果是否正确
bool verifyResult(float *C1, float *C2, int n) {float epsilon = 1e-4f;for (int i = 0; i < n * n; i++) {if (fabs(C1[i] - C2[i]) > epsilon) {printf("验证失败! 索引 %d: CPU=%f, GPU=%f\n", i, C1[i], C2[i]);return false;}}return true;
}int main() {// 分配主机内存size_t size = N * N * sizeof(float);float *h_A = (float*)malloc(size);float *h_B = (float*)malloc(size);float *h_C_cpu = (float*)malloc(size);float *h_C_gpu_simple = (float*)malloc(size);float *h_C_gpu_shared = (float*)malloc(size);// 初始化矩阵initMatrix(h_A, N, 1);initMatrix(h_B, N, 2);// 分配设备内存float *d_A, *d_B, *d_C;cudaMalloc(&d_A, size);cudaMalloc(&d_B, size);cudaMalloc(&d_C, size);// 拷贝数据到设备cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);// 设置CUDA网格和块维度dim3 threadsPerBlock(BLOCK_SIZE, BLOCK_SIZE);dim3 blocksPerGrid((N + BLOCK_SIZE - 1) / BLOCK_SIZE, (N + BLOCK_SIZE - 1) / BLOCK_SIZE);// 计时变量clock_t start, end;cudaEvent_t gpuStart, gpuStop;cudaEventCreate(&gpuStart);cudaEventCreate(&gpuStop);float cpuTime, gpuSimpleTime, gpuSharedTime;// 运行CPU版本并计时printf("运行CPU矩阵乘法...\n");start = clock();matrixMulCPU(h_A, h_B, h_C_cpu, N);end = clock();cpuTime = ((float)(end - start)) / CLOCKS_PER_SEC;printf("CPU时间: %.3f 秒\n", cpuTime);// 运行简单GPU版本并计时printf("运行简单GPU矩阵乘法...\n");cudaEventRecord(gpuStart);matrixMulSimple<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);cudaEventRecord(gpuStop);cudaEventSynchronize(gpuStop);cudaMemcpy(h_C_gpu_simple, d_C, size, cudaMemcpyDeviceToHost);cudaEventElapsedTime(&gpuSimpleTime, gpuStart, gpuStop);gpuSimpleTime /= 1000.0f; // 转换为秒printf("简单GPU时间: %.3f 秒, 加速比: %.2fx\n", gpuSimpleTime, cpuTime / gpuSimpleTime);// 验证简单GPU版本结果printf("验证简单GPU结果...\n");if (verifyResult(h_C_cpu, h_C_gpu_simple, N)) {printf("简单GPU结果正确!\n");} else {printf("简单GPU结果错误!\n");}// 运行共享内存GPU版本并计时printf("运行共享内存GPU矩阵乘法...\n");cudaEventRecord(gpuStart);matrixMulShared<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);cudaEventRecord(gpuStop);cudaEventSynchronize(gpuStop);cudaMemcpy(h_C_gpu_shared, d_C, size, cudaMemcpyDeviceToHost);cudaEventElapsedTime(&gpuSharedTime, gpuStart, gpuStop);gpuSharedTime /= 1000.0f; // 转换为秒printf("共享内存GPU时间: %.3f 秒, 加速比: %.2fx\n", gpuSharedTime, cpuTime / gpuSharedTime);// 验证共享内存GPU版本结果printf("验证共享内存GPU结果...\n");if (verifyResult(h_C_cpu, h_C_gpu_shared, N)) {printf("共享内存GPU结果正确!\n");} else {printf("共享内存GPU结果错误!\n");}// 清理资源free(h_A);free(h_B);free(h_C_cpu);free(h_C_gpu_simple);free(h_C_gpu_shared);cudaFree(d_A);cudaFree(d_B);cudaFree(d_C);cudaEventDestroy(gpuStart);cudaEventDestroy(gpuStop);return 0;
}

如何理解 cuda 代码的索引计算？

    int row = blockIdx.y * BLOCK_SIZE + threadIdx.y;int col = blockIdx.x * BLOCK_SIZE + threadIdx.x;

在核函数内部，CUDA提供了几个关键的内置变量来标识当前线程的位置：

threadIdx.x, threadIdx.y, threadIdx.z: 当前线程在其所属线程块中的索引。

blockIdx.x, blockIdx.y, blockIdx.z: 当前线程块在整个网格中的索引。

blockDim.x, blockDim.y, blockDim.z: 线程块的维度（即，每个块有多少个线程，你在<<<>>>中指定的threadsPerBlock）。

为什么共享内存可以加速矩阵乘法？

共享内存（Shared Memory）之所以能极大提升性能，主要是因为它解决了GPU计算中最主要的性能瓶颈：全局内存（Global Memory）的访问延迟和带宽限制。

我们可以通过一个生动的比喻和深入的技术分析来理解这一点。

1. 核心思想：数据复用 (Data Reuse)

许多算法，尤其是像矩阵乘法、卷积这样的科学计算，存在大量的数据复用。即同一份数据需要被使用多次。

在矩阵乘法中：矩阵A的一行需要与矩阵B的每一列进行点乘。同样，矩阵B的一列也需要与矩阵A的每一行进行点乘。

2. 没有共享内存的问题：“直接”访问全局内存

如果使用简单的核函数（如 matrixMulSimple），计算过程是这样的：

每个线程读取矩阵A的一整行和矩阵B的一整列。
每个数据元素都是从全局内存中读取的。
对于 N x N 的矩阵，每个线程需要读取 2N 次全局内存，总共是 2N³ 次全局内存访问。

全局内存的特点：

速度慢：延迟非常高（需要几百甚至上千个时钟周期）。
带宽高但宝贵：虽然总带宽很大，但被成千上万个线程争抢，容易成为瓶颈。
耗能高：访问全局内存是GPU中最耗能的操作之一。

这导致了巨大的浪费：例如，矩阵A的第i行会被计算C[i][j]（j从0到N-1）的所有线程重复读取N次。每次读取都要经历漫长且耗能的全局内存访问。

3. 共享内存的解决方案：充当高速缓存 (On-Chip Cache)

共享内存是一块位于GPU芯片上（On-Chip）的高速、低延迟的存储器，由同一个线程块（Block）内的所有线程共享。

优化后的流程（如 matrixMulShared）：

协作加载 (Cooperative Loading)：一个线程块内的所有线程协作将计算所需的一小块数据（Tile）从慢速的全局内存加载到快速的共享内存中。
- 加载次数：(N / BLOCK_SIZE) * 2 * (BLOCK_SIZE²) 次。
高速访问：线程块内的所有计算都直接从共享内存中读取所需数据。
- 共享内存的延迟比全局内存低**_{100倍**，带宽高**}10倍**。
同步：使用 __syncthreads() 确保所有线程都完成数据加载后，才开始计算，防止数据竞争。

这就好比：

没有共享内存：你要做一顿大餐，每次需要一种食材，你就开车去很远的大型超市（全局内存）买回来，用完再去买下一种。大部分时间都花在路上了。
使用共享内存：你先和你的团队（线程块）一起开车去超市，把做这几道菜需要的所有食材（数据块）一次性大批量采购回来，放在厨房的台面上（共享内存）。接下来整个烹饪过程（计算）都在台面上快速完成，再也无需往返超市。

4. 性能提升的关键总结

大幅减少全局内存访问次数：
- 简单版本：O(N³) 次访问。
- 共享内存版本：O(N²) 次访问（每个数据元素只从全局内存加载一次）。
- 访问次数降低了一个数量级（从N³到N²），这是性能提升的最主要原因。
利用共享内存的高速度：
- 后续成百上千次的数据访问都发生在纳秒级别的共享内存上，而不是微秒级别的全局内存上。
促进内存访问合并 (Memory Coalescing)：
- 线程协作加载时，可以将对全局内存的访问模式从“分散”的变为“连续”的。连续的内存访问可以被合并成一次大的事务，从而更高效地利用全局内存的高带宽。

一个简单的计算对比

假设 N=1024, BLOCK_SIZE=32：

简单核函数：
- 总全局内存访问量 ≈ 2 * 1024³ = ~21.5亿次访问。
共享内存核函数：
- 每个数据元素（A和B中的每个值）只被从全局内存加载一次到共享内存。
- 总全局内存访问量 ≈ 2 * 1024² = ~210万次访问。
- 然后，在共享内存中会发生 1024³ 次访问，但这些访问速度极快。

可以看到，对全局内存的访问次数减少了 1024倍！尽管增加了同步和共享内存管理的开销，但净收益是巨大的。

注意事项

共享内存并非“银弹”，使用它需要遵循一些规则：

Bank Conflict：共享内存被组织成多个bank。如果同一个warp内的多个线程同时访问同一个bank的不同地址，就会发生bank conflict，导致访问被序列化，降低性能。优秀的CUDA程序员会精心设计数据在共享内存中的布局和访问模式来避免这一点。
同步开销：不必要或过多的 __syncthreads() 会影响性能。
适用范围：最适合存在显著数据复用的场景。如果每个数据只使用一次，那么使用共享内存反而会增加额外开销，降低性能。

总而言之，共享内存优化通过让线程块协作地将数据从慢速的全局内存“缓存”到快速的片上内存，并最大限度地复用这些数据，从而避免了重复、昂贵的长距离内存访问，这是其性能提升的根本原因。这是CUDA编程中最重要、最经典的优化手段之一。