CUDA 矩阵乘法 优化版本

//矩阵乘法的优化版本

#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;




}