//矩阵乘法的优化版本
#include<cuda_runtime.h>
#include<device_launch_parameters.h>
#include<stdlib.h>
#include<stdio.h>
#include<assert.h>
#include<math.h>
#define SHMEM_SIZE 16*16*4 //为方便起见,静态shmem计算 16x16矩阵
__global__ void tiledMatrixMul(int*a, int*b, int*c, int n,int tile_size) {
__shared__ int A[SHMEM_SIZE];
__shared__ int B[SHMEM_SIZE];//两个静态的共享内存
int tx = threadIdx.x;
int ty = threadIdx.y;
int bx = blockIdx.x;
int by = blockIdx.y;//缩短参数
int row = by * tile_size + ty;
int col = bx * tile_size + tx;//计算该线程的全局行和列位置
int temp_val = 0;//被写入元素的中间和
//将图块扫过整个矩阵
for (int i = 0; i < (n / tile_size); i++) {
/*
线程块中每个线程都将一个元素加载到共享内存中
共享内存中的元素位置对应于线程的线程块中的位置
(比如,线程[0,0]加载为A[0*tile_size+0]和B[0*tile_size+0]
索引参数说明
A:row*n:索引此线程的全局行(循环不变)
i*tile_size:每次迭代索引新的一组列
tx:索引该集中的列
B:i*tile_size*n:每次迭代索引下一组行
ty*n:索引该集合中的行
col:索引全局列(循环不变)
*/
A[(ty*tile_size) + tx] = a[row*n + (i*tile_size + tx)];
B[(ty*tile_size) + tx] = b[(i*tile_size*n + ty * n) + col];
__syncthreads();//同步 确保所有线程都已加载数据,然后再继续
//计算此图块的所有临时值
for (int j = 0; j < tile_size; j++) {
temp_val += A[(ty*tile_size) + j] * B[(j*tile_size) + tx];
}
__syncthreads();//同步 确保某些线程不会继续运行并占用当前共享内存值
}
c[(row*n) + col] = temp_val;
}
void check_answer(int*a, int*b, int*c, int n) {
int tmp;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
tmp = 0;
for (int k = 0; k < n; k++) {
tmp += a[i*n + k] * b[k*n + j];
}
assert(tmp == c[i*n + j]);
}
}
}
void init_matrix(int *a, int n) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
a[i*n + j] = rand() % 10;
}
}
}
int main() {
int n = 1 << 10;
size_t bytes = n * n*sizeof(int);
int *h_a, *h_b, *h_c;
int *d_a, *d_b, *d_c;
h_a = (int*)malloc(bytes);
h_b = (int*)malloc(bytes);
h_c = (int*)malloc(bytes);
cudaMalloc(&d_a, bytes);
cudaMalloc(&d_b, bytes);
cudaMalloc(&d_c, bytes);
init_matrix(h_a, n);
init_matrix(h_b, n);
cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice);
int BLOCK_SIZE = 16;
int GRID_SIZE = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
dim3 grid(GRID_SIZE, GRID_SIZE);
dim3 threads(BLOCK_SIZE, BLOCK_SIZE);
tiledMatrixMul << <grid, threads >> > (d_a, d_b, d_c, n, BLOCK_SIZE);
cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);
check_answer(h_a,h_b, h_c, n);
free(h_a);
free(h_b);
free(h_c);
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);
printf("completed successfully\n");
return 0;
}