cuda学习（三）矩阵乘法
一、方法一
#include <iostream>
#include <cuda_runtime.h>
__global__ void matrixMultiply(int *A, int *B, int *C, int N) {
// 计算每个线程负责的矩阵元素位置
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < N && col < N) {
int value = 0;
// 计算矩阵 C 的元素 C[row][col]
for (int k = 0; k < N; ++k) {
value += A[row * N + k] * B[k * N + col];
}
C[row * N + col] = value;
}
}
int main() {
int N = 1024; // 矩阵大小 N x N
int size = N * N * sizeof(int); // 每个矩阵的字节大小
// 为矩阵 A, B, C 在主机上分配内存
int *h_A = new int[N * N];
int *h_B = new int[N * N];
int *h_C = new int[N * N];
// 初始化矩阵 A 和 B
for (int i = 0; i < N * N; ++i) {
h_A[i] = 1; // 矩阵 A 初始化为 1
h_B[i] = 2; // 矩阵 B 初始化为 2
}
// 为设备上的矩阵 A, B, C 分配内存
int *d_A, *d_B, *d_C;
cudaMalloc(&d_A, size);
cudaMalloc(&d_B, size);
cudaMalloc(&d_C, size);
// 将矩阵 A 和 B 从主机复制到设备
cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);
// 设置线程和块的维度
dim3 threadsPerBlock(16, 16); // 每个块 16x16 的线程
dim3 blocksPerGrid((N + 15) / 16, (N + 15) / 16); // 网格的块的数量
// 启动内核进行矩阵乘法
matrixMultiply<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);
// 检查内核启动是否有错误
cudaError_t error = cudaGetLastError();
if (error != cudaSuccess) {
std::cerr << "CUDA kernel launch failed: " << cudaGetErrorString(error) << std::endl;
return -1;
}
// 等待设备上的所有计算完成
cudaDeviceSynchronize();
// 将结果从设备复制回主机
cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);
// 打印部分结果验证
for (int i = 0; i < 10; ++i) { // 打印前 10 个结果
std::cout << "C[" << i << "] = " << h_C[i] << std::endl;
}
// 释放设备内存
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);
// 释放主机内存
delete[] h_A;
delete[] h_B;
delete[] h_C;
return 0;
}
方法一和矩阵加法类似，不再赘述
二、使用共享内存来优化矩阵乘法
#include <iostream>
#include <cuda_runtime.h>
#define TILE_SIZE 16 // 线程块的大小，16x16 的块
__global__ void matrixMultiplyWithSharedMemory(int *A, int *B, int *C, int N) {
__shared__ int s_A[TILE_SIZE][TILE_SIZE]; // 用于存储矩阵 A 的共享内存
__shared__ int s_B[TILE_SIZE][TILE_SIZE]; // 用于存储矩阵 B 的共享内存
int row = blockIdx.y * TILE_SIZE + threadIdx.y; // 当前线程处理的矩阵行
int col = blockIdx.x * TILE_SIZE + threadIdx.x; // 当前线程处理的矩阵列
int value = 0;
for (int i = 0; i < (N / TILE_SIZE); ++i) {
// 将 A 和 B 的子块加载到共享内存
s_A[threadIdx.y][threadIdx.x] = A[row * N + (i * TILE_SIZE + threadIdx.x)];
s_B[threadIdx.y][threadIdx.x] = B[(i * TILE_SIZE + threadIdx.y) * N + col];
__syncthreads(); // 确保所有线程都完成了共享内存加载
// 执行子块的乘法
for (int j = 0; j < TILE_SIZE; ++j) {
value += s_A[threadIdx.y][j] * s_B[j][threadIdx.x];
}
__syncthreads(); // 等待所有线程完成本次子块计算
}
if (row < N && col < N) {
C[row * N + col] = value; // 将计算结果存储到矩阵 C 中
}
}
int main() {
int N = 1024; // 矩阵的大小 N x N
int size = N * N * sizeof(int); // 每个矩阵的字节大小
// 在主机上分配矩阵 A, B, C 的内存
int *h_A = new int[N * N];
int *h_B = new int[N * N];
int *h_C = new int[N * N];
// 初始化矩阵 A 和 B
for (int i = 0; i < N * N; ++i) {
h_A[i] = 1; // 矩阵 A 初始化为 1
h_B[i] = 2; // 矩阵 B 初始化为 2
}
// 为设备上的矩阵 A, B, C 分配内存
int *d_A, *d_B, *d_C;
cudaMalloc(&d_A, size);
cudaMalloc(&d_B, size);
cudaMalloc(&d_C, size);
// 将矩阵 A 和 B 从主机复制到设备
cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);
// 设置线程块和网格的维度
dim3 threadsPerBlock(TILE_SIZE, TILE_SIZE); // 每个线程块 16x16 的线程
dim3 blocksPerGrid((N + TILE_SIZE - 1) / TILE_SIZE, (N + TILE_SIZE - 1) / TILE_SIZE); // 网格中块的数量
// 启动内核进行矩阵乘法
matrixMultiplyWithSharedMemory<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);
// 检查内核启动是否有错误
cudaError_t error = cudaGetLastError();
if (error != cudaSuccess) {
std::cerr << "CUDA kernel launch failed: " << cudaGetErrorString(error) << std::endl;
return -1;
}
// 等待设备上的所有计算完成
cudaDeviceSynchronize();
// 将结果从设备复制回主机
cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);
// 打印部分结果验证
for (int i = 0; i < 10; ++i) { // 打印前 10 个结果
std::cout << "C[" << i << "] = " << h_C[i] << std::endl;
}
// 释放设备内存
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);
// 释放主机内存
delete[] h_A;
delete[] h_B;
delete[] h_C;
return 0;
}
先贴上代码。
第一步：申请主机内存
int N = 1024; // 矩阵的大小 N x N
int size = N * N * sizeof(int); // 每个矩阵的字节大小
// 在主机上分配矩阵 A, B, C 的内存
int *h_A = new int[N * N];
int *h_B = new int[N * N];
int *h_C = new int[N * N];
第二步：设备内存分配
cudaMalloc(&d_A, size); // 在设备（GPU）上为矩阵 A 分配内存
cudaMalloc(&d_B, size); // 在设备（GPU）上为矩阵 B 分配内存
cudaMalloc(&d_C, size); // 在设备（GPU）上为矩阵 C 分配内存
使用 cudaMalloc 在 GPU 上为矩阵 A、B 和 C 分配内存。每个矩阵的大小为 N * N * sizeof(int) 字节。
第三步：初始化数据
for (int i = 0; i < N * N; ++i) {
h_A[i] = 1; // 矩阵 A 初始化为 1
h_B[i] = 2; // 矩阵 B 初始化为 2
}
这部分代码初始化了矩阵 A 和矩阵 B。A 的每个元素都设置为 1，B 的每个元素都设置为 2。这里为了简单起见，矩阵值的初始化较为简单。
第四步：将数据从主机复制到设备
cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);
使用 cudaMemcpy 将主机内存中的矩阵 A 和 B 复制到 GPU 的设备内存中。这样，GPU 才能够访问这些数据进行计算。
第五步：内核函数编写
__global__ void matrixMultiplyWithSharedMemory(int *A, int *B, int *C, int N) {
__shared__ int s_A[TILE_SIZE][TILE_SIZE]; // 用于存储矩阵 A 的共享内存
__shared__ int s_B[TILE_SIZE][TILE_SIZE]; // 用于存储矩阵 B 的共享内存
int row = blockIdx.y * TILE_SIZE + threadIdx.y; // 当前线程处理的矩阵行
int col = blockIdx.x * TILE_SIZE + threadIdx.x; // 当前线程处理的矩阵列
int value = 0;
for (int i = 0; i < (N / TILE_SIZE); ++i) {
// 将 A 和 B 的子块加载到共享内存
s_A[threadIdx.y][threadIdx.x] = A[row * N + (i * TILE_SIZE + threadIdx.x)];
s_B[threadIdx.y][threadIdx.x] = B[(i * TILE_SIZE + threadIdx.y) * N + col];
__syncthreads(); // 确保所有线程都完成了共享内存加载
// 执行子块的乘法
for (int j = 0; j < TILE_SIZE; ++j) {
value += s_A[threadIdx.y][j] * s_B[j][threadIdx.x];
}
__syncthreads(); // 等待所有线程完成本次子块计算
}
if (row < N && col < N) {
C[row * N + col] = value; // 将计算结果存储到矩阵 C 中
}
}
__shared__ 关键字表示该变量是在 共享内存 中分配的。共享内存是每个线程块（block）内部的高速内存，线程块内的所有线程都能访问共享内存。相比于全局内存（global memory），共享内存的访问速度快很多，因此在计算中合理利用共享内存可以大大提升性能。这里使用了 2 个 TILE_SIZE x TILE_SIZE 的二维数组 s_A 和 s_B 来存储矩阵 A 和 B 的子块。这是因为 GPU 上的线程块（block）会负责处理矩阵的一个子块。
blockIdx.x 和 blockIdx.y 是当前线程块的网格坐标，表示线程块在网格中的位置。threadIdx.x 和 threadIdx.y 是当前线程在所在线程块中的位置。TILE_SIZE 是线程块的大小，这里设置为 16，表示每个线程块是 16x16 的线程矩阵。row 和 col 用来计算每个线程对应矩阵 C 中的元素位置。每个线程负责计算结果矩阵中的一个元素。
A[row * N + (i * TILE_SIZE + threadIdx.x)] 计算当前线程要从矩阵 A 中读取的元素的位置。B[(i * TILE_SIZE + threadIdx.y) * N + col] 计算当前线程要从矩阵 B 中读取的元素的位置。
第六步：内核调用与执行
dim3 threadsPerBlock(TILE_SIZE, TILE_SIZE); // 每个线程块 16x16 的线程
dim3 blocksPerGrid((N + TILE_SIZE - 1) / TILE_SIZE, (N + TILE_SIZE - 1) / TILE_SIZE); // 网格中块的数量
matrixMultiplyWithSharedMemory<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);
threadsPerBlock 设置每个线程块的线程数为 TILE_SIZE x TILE_SIZE，即每个线程块包含 16x16 个线程，共 256 个线程。blocksPerGrid 根据矩阵的大小和线程块的大小计算出所需的线程块数量，保证每个矩阵元素都有一个线程负责计算。内核函数 matrixMultiplyWithSharedMemory 被启动并执行。
后面重复过程不再赘述。
三、一些问题的思考
1.在使用共享内存时，线程块内的所有线程都访问共享内存。如果多个线程同时访问共享内存的相同位置，会不会造成冲突？如果会，应该如何避免？
多个线程同时访问共享内存的同一位置会导致冲突，称为 共享内存银行冲突。为了避免这种情况，CUDA 的共享内存被分为多个 银行，每个银行可以并行访问。但如果多个线程访问同一个银行的同一位置，就会发生冲突，影响性能。
要避免共享内存银行冲突，可以通过调整内存访问模式，使得每个线程访问不同的内存位置。例如，通过改变 TILE_SIZE 的大小或者通过更精细的内存布局，确保不同线程访问不同的银行。
2.如果矩阵的维度 N 不是 TILE_SIZE 的倍数，如何处理那些无法被 TILE_SIZE 整除的部分？当前代码是否能有效处理这种情况？
可以修改内核代码，在加载子块数据时，确保矩阵的边界被正确处理。例如，在加载 A 和 B 子块时，检查 row < N 和 col < N，并在矩阵 C 中写回计算结果时，确保只更新有效元素。
3.矩阵 A 和 B 的加载顺序是否最优？是否有可能优化内存访问模式，减少不必要的内存访问延迟？
对于矩阵 B，可以考虑对其进行 转置，使得每个线程加载相邻的内存位置，从而提高内存访问效率。转置后的 B 会使得线程块在加载数据时，线程按行访问，这可以减少不必要的内存访问延迟。对于矩阵 A，可以保持当前加载顺序，但对于矩阵 B，转置可能会提升性能。
4.内核的 for (int i = 0; i < (N / TILE_SIZE); ++i) 语句假设 N 是 TILE_SIZE 的倍数。若 N 不能整除 TILE_SIZE，代码是否会遇到问题？如何处理这个问题？
为了处理非 TILE_SIZE 倍数的情况，可以计算子块数量时向上取整（ceil）。例如，可以使用 (N + TILE_SIZE - 1) / TILE_SIZE 计算需要的子块数量，确保覆盖矩阵的所有部分。
5.当前代码每个线程块都需要从全局内存加载矩阵 A 和 B 的子块到共享内存，是否有可能减少内存传输开销？
减少全局内存访问的次数，可以通过将矩阵 A 和 B 中的多个子块存储在共享内存中并重用，减少不必要的内存传输。另外，尝试通过合理设计 内存对齐 和 内存访问模式 来进一步减少内存传输的开销。
5.可不可以将整个A和B矩阵整体传入共享内存然后每次取一个方块？
CUDA 中的共享内存是有限制的。对于每个线程块（block），共享内存的大小通常是每个块的 16KB 或者 48KB，具体大小取决于 GPU 架构和计算能力（Compute Capability）。如果矩阵的大小 N 很大（例如 1024 x 1024），即使是存储一个完整的矩阵也会占用较多的共享内存。例如，N = 1024 时，每个矩阵 A 和 B 的内存大小为 1024 x 1024 x sizeof(int)，即 1024 x 1024 x 4 = 4MB，这显然远远超过了单个线程块的共享内存容量。因此，直接将整个矩阵加载到共享内存并不是一个可行的选择，除非矩阵的尺寸非常小，或者使用多个线程块来分担内存工作。
更多推荐
17
0- 0
已为社区贡献1条内容
所有评论(0)