CUDA矩阵乘法

#include "stdio.h"
#include<iostream>
#include <cuda.h>
#include <cuda_runtime.h>
#include <math.h>

#define TILE_SIZE 2

//Matrix multiplication using non shared kernel
__global__ void gpu_Matrix_Mul_nonshared(float *d_a, float *d_b, float *d_c, const int size)
{
	int row, col;
	col = TILE_SIZE * blockIdx.x + threadIdx.x;
	row = TILE_SIZE * blockIdx.y + threadIdx.y;

	for (int k = 0; k < size; k++)
	{
		d_c[row*size + col] += d_a[row * size + k] * d_b[k * size + col];
	}
}

// Matrix multiplication using shared kernel
__global__ void gpu_Matrix_Mul_shared(float *d_a, float *d_b, float *d_c, const int size)
{
	int row, col;
	//Defining Shared Memory
	__shared__ float shared_a[TILE_SIZE][TILE_SIZE];
	__shared__ float shared_b[TILE_SIZE][TILE_SIZE];
	col = TILE_SIZE * blockIdx.x + threadIdx.x;
	row = TILE_SIZE * blockIdx.y + threadIdx.y;

	for (int i = 0; i < size / TILE_SIZE; i++) 
	{
		shared_a[threadIdx.y][threadIdx.x] = d_a[row* size + (i*TILE_SIZE + threadIdx.x)];
		shared_b[threadIdx.y][threadIdx.x] = d_b[(i*TILE_SIZE + threadIdx.y) * size + col];
		__syncthreads(); 
		for (int j = 0; j<TILE_SIZE; j++)
			d_c[row*size + col] += shared_a[threadIdx.y][j] * shared_b[j][threadIdx.x];
		__syncthreads(); 
	}
}

int main()
{
	const int size = 4;
	//Define Host Array
	float h_a[size][size], h_b[size][size],h_result[size][size];
	//Defining device Array
	float *d_a, *d_b, *d_result; 
	//Initialize host Array
	for (int i = 0; i<size; i++)
	{
		for (int j = 0; j<size; j++)
		{
			h_a[i][j] = i;
			h_b[i][j] = j;
		}
	}

	cudaMalloc((void **)&d_a, size*size*sizeof(int));
	cudaMalloc((void **)&d_b, size*size * sizeof(int));
	cudaMalloc((void **)&d_result, size*size* sizeof(int));

	//copy host array to device array
	cudaMemcpy(d_a, h_a, size*size* sizeof(int), cudaMemcpyHostToDevice);
	cudaMemcpy(d_b, h_b, size*size* sizeof(int), cudaMemcpyHostToDevice);
	
	//Define grid and block dimensions
	dim3 dimGrid(size / TILE_SIZE, size / TILE_SIZE, 1);
	dim3 dimBlock(TILE_SIZE, TILE_SIZE, 1);
	//gpu_Matrix_Mul_nonshared << <dimGrid, dimBlock >> > (d_a, d_b, d_result, size);
	gpu_Matrix_Mul_shared << <dimGrid, dimBlock >> > (d_a, d_b, d_result, size);

	cudaMemcpy(h_result, d_result, size*size * sizeof(int),	cudaMemcpyDeviceToHost);
	printf("The result of Matrix multiplication is: \n");
	
	for (int i = 0; i< size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			printf("%f   ", h_result[i][j]);
		}
		printf("\n");
	}

	cudaFree(d_a);
	cudaFree(d_b);
	cudaFree(d_result);
	return 0;
}

矩阵乘法中同样的数据能被使用多次，所以这是共享内存的一个理想用例。在本小节中，我们写出两个独立的内核，它们将会分别使用和不使用共享内存。你可以通过比较这两个内核的执行(速度)，来了解共享内存对程序的性能提升效果。我们先开始写一个不使用共享内存的内核，每个块中的线程将计算这些小方块中（特定结果矩阵)的元素。进行矩阵乘法的总块数量将是用矩阵的大小除以小方块的大小(TILE_SIZE)计算得到。
如果你能理解切分问题，则计算(每个线程对应)结果矩阵中的行列索引值就非常容易了。它很像我们之前计算(全局索引)的方式，只是这里的 blockDim.x和 blockDim.y等于TILE_SIZE的大小。现在，结果矩阵中的每个元素都将是第一个矩阵的(对应)行和第二个矩阵的(对应)列进行向量点乘的结果。因为两个矩阵都是变量size * size这么大的方阵，所以需要是size个元素这么大的两个向量进行点乘。所以(每个线程)通过内核函数代码里的for循环从О循环到size，进行了多次乘加。
为了能够在两个矩阵中一个一个的索引元素(计算元素的索引)，可以将矩阵看成是在存储器中按照行主序的方式线性存储的。行主序的意思是，第一行的元素一个接一个地连续在存储器中排列，然后下一行再继续(这样排列)在前一行的后面。如下所示:

每个元素的线性索引可以这样计算:用它的行号乘以矩阵的宽度，再加上它的列号即可。所以，对于M,o来说，它的线性索引值是2。因为它的行号是1，矩阵的宽度是2，它的列号是0。这种方法用来计算两个源矩阵中元素的索引。
为了计算结果矩阵中在[row,col]位置的元素的值，需要前一个矩阵中的一行和后一个矩阵的一列进行点乘。内核里分别用了rowsize+k的下标和ksize+col的下标来索引这一行和一列中的元素。
要计算最终的目标矩阵中的[row,col]的元素，需要将第一个矩阵中的第row行和第二个矩阵中的第col列进行向量点乘。而这向量点乘则需要循环多组对应的点，分别对来自第一个矩阵从2D索引线性化到1D索引后的坐标为rowsize+k的点，和来自第二个矩阵的同样1D化坐标为ksize+col处的点进行乘法累加操作，其中k是每次循环时候的变量。这样就完成了向量点乘,以及最后的矩阵乘法。
用共享内存中的大小为TILE_SIZE * TILE_SIZE的两个方块来保存需要重用的数据。就像你之前看到的那样的方式来计算行和列索引。首次for循环的时候，先填充共享内存，然后调用__syncthreads()，从而确保后续从共享内存的读取只会发生在块中的所有线程都完成填充共享内存之后。然后循环的后面部分则是计算(部分）向量点乘。因为(大量的重复读取）都只发生在共享内存上，从而显著地降低对全局内存的访存流量，进一步地能提升大矩阵维度下的乘法程序的性能。
当定义了主机和设备指针变量并为它们分配空间后，主机上的数组被填充了一些任意的数字(注意不是随机，因为是写死的)，然后这些数组被复制到显存上，这样它们才能作为参数传输给内核函数。接着使用dim3结构体定义Grid中的块和块中的线程形状，具体的形状信息是之前几行算好的，你可以调用这两个内核中的任何一个(使用共享内存的和不使用共享内存的)，然后将结果复制回主机来。