一个用GPU写的矩阵乘法

这个过程相当繁琐,个人认为有优化的可能:

先说一下思路,矩阵相乘A矩阵乘B矩阵相当于A矩阵和B矩阵的转置做内积.所以我就先把B矩阵做了转置,再做内积.其中有两个核函数是在主函数中执行的,先执行转置,再执行乘法.再乘法函数中又嵌套了一个内积函数.这样充分的利用了并行化.

如图所示:

以3*3矩阵为例.我先开3*3个线程做内积运算,然后在每个线程中又开了1*3个线程做内积运算.其中求和部分我没有用并行方式求和而是使用串行方式求和的,个人认为这是一个可以优化的点.还有一个优化点就是,个人认为不必要算矩阵的转置就可以实现内积,只不过在乘法中算法变一下就行了.以后有机会讨论优化问题.这样我就开了3*3*3等于27个线程.这样求出的矩阵乘法.

一会看看其他人的算法思路是否有借鉴之处,这是本人的GPU求矩阵的乘法的思路.思想很简单,但是不是很成熟,有更简单的算法一起讨论学习.

代码使用:

输入一个数为矩阵的维度,代码会自动生成两个n*n的随机矩阵.然后相乘得到另外的矩阵.

实验结果:

请输入矩阵的维度:
3

CPU内存数据h_a:
h_a[0][0] = 2	h_a[0][1] = 4	h_a[0][2] = 10	
h_a[1][0] = 9	h_a[1][1] = 7	h_a[1][2] = 3	
h_a[2][0] = 8	h_a[2][1] = 10	h_a[2][2] = 9	
CPU内存数据h_b:
h_b[0][0] = 2	h_b[0][1] = 0	h_b[0][2] = 8	
h_b[1][0] = 7	h_b[1][1] = 4	h_b[1][2] = 8	
h_b[2][0] = 6	h_b[2][1] = 1	h_b[2][2] = 6	
CPU内存数据h_bt:
h_bt[0][0] = 2	h_bt[0][1] = 7	h_bt[0][2] = 6	
h_bt[1][0] = 0	h_bt[1][1] = 4	h_bt[1][2] = 1	
h_bt[2][0] = 8	h_bt[2][1] = 8	h_bt[2][2] = 6	
GPU内存数据:
h_c[0][0] = 92	h_c[0][1] = 26	h_c[0][2] = 108	
h_c[1][0] = 85	h_c[1][1] = 31	h_c[1][2] = 146	
h_c[2][0] = 140	h_c[2][1] = 49	h_c[2][2] = 198	
运行结束

 代码:

#include <cuda_runtime.h>
#include <iostream>
#include <stdlib.h>
#include <time.h>

__global__ void zhuanshi(int *d_b,int *d_bt,int n)
{
	int ITdx = threadIdx.x;
	int IBdx = blockIdx.x;

	d_bt[ITdx * n + IBdx] = d_b[IBdx * n + ITdx];
}
__global__ void neiji(int *d_a,int *d_bt,int *d_c,int *d_data,int ICTdx,int ICBdx,int n)
{
	int INTdx = threadIdx.x;
	int temp = 0;

	d_data[(ICTdx * n + ICBdx) * n + INTdx] = d_a[ICTdx * n + INTdx] * d_bt[ ICBdx * n + INTdx];

	__syncthreads();
	for(int i = 0; i < n; ++i)
	{
		temp += d_data[(ICTdx * n + ICBdx) * n + i];
	}
	d_c[ICTdx * n + ICBdx] = temp;

}
__global__ void chengfa(int *d_a,int *d_bt,int *d_c,int *d_data,int n)
{

	int ICTdx = threadIdx.x;
	int ICBdx = blockIdx.x;
	neiji<<<1,n>>>(d_a,d_bt,d_c,d_data,ICTdx,ICBdx,n);
	__syncthreads();

}

int main()
{
	int blag = 1;//标志位
	int n = 0;
	/******判断输入数据是否合法************/
	do{
		std::cout << "请输入矩阵的维度:" << std::endl;
		std::cin >> n;
		if(n <= 0)
		{
			std::cout << "你输入的矩阵维度有误,请重新输入!" << std::endl;
		}else{
			blag = 0;
		}
	}while(blag);

	/*******申请主机内存*********/
	int *h_a = (int*)malloc(sizeof(int) * n * n);
	int *h_b = (int*)malloc(sizeof(int) * n * n);
	int *h_c = (int*)malloc(sizeof(int) * n * n);
	int *h_bt = (int*)malloc(sizeof(int) * n * n);
	/*******初始化主机内存数据********/
	srand(time(NULL));//设置随机数值
	for(int i = 0; i < n * n; ++i)
	{
		h_a[i] = rand() % 11;
		h_b[i] = rand() % 11;
		h_c[i] = 0;
		h_bt[i] = 0;
	}

	/*******申请设备内存*******/
	int *d_a,*d_b,*d_c,*d_bt,*d_data;
	cudaMalloc((void**)&d_a,sizeof(int) * n * n);
	cudaMalloc((void**)&d_b,sizeof(int) * n * n);
	cudaMalloc((void**)&d_c,sizeof(int) * n * n);
	cudaMalloc((void**)&d_bt,sizeof(int) * n * n);
	cudaMalloc((void**)&d_data,sizeof(int) * n * n * n);

	/******主机内存数据复制到设备内存中************/
	cudaMemcpy(d_a,h_a,sizeof(int) * n * n,cudaMemcpyHostToDevice);
	cudaMemcpy(d_b,h_b,sizeof(int) * n * n,cudaMemcpyHostToDevice);

	/*******执行核函数******/
	zhuanshi<<<n,n>>>(d_b,d_bt,n);
	chengfa<<<n,n>>>(d_a,d_bt,d_c,d_data,n);

	/*****设备内存数据复制到主机内存中*****/
	cudaMemcpy(h_bt,d_bt,sizeof(int) * n * n,cudaMemcpyDeviceToHost);
	cudaMemcpy(h_c,d_c,sizeof(int) * n * n,cudaMemcpyDeviceToHost);

	std::cout << "CPU内存数据h_a:" << std::endl;
	for(int i = 0; i < n; ++i)
	{
		for(int j = 0; j < n; ++j)
		{
			printf("h_a[%d][%d] = %d\t",i,j,h_a[n * i + j]);
		}
		printf("\n");
	}
	std::cout << "CPU内存数据h_b:" << std::endl;
	for(int i = 0; i < n; ++i)
	{
		for(int j = 0; j < n; ++j)
		{
			printf("h_b[%d][%d] = %d\t",i,j,h_b[n * i + j]);
		}
		printf("\n");
	}
	std::cout << "CPU内存数据h_bt:" << std::endl;
	for(int i = 0; i < n; ++i)
	{
		for(int j = 0; j < n; ++j)
		{
			printf("h_bt[%d][%d] = %d\t",i,j,h_bt[n * i + j]);
		}
		printf("\n");
	}
	std::cout << "GPU内存数据:" << std::endl;
	for(int i = 0; i < n; ++i)
	{
		for(int j = 0; j < n; ++j)
		{
			printf("h_c[%d][%d] = %d\t",i,j,h_c[n * i + j]);
		}
		printf("\n");
	}

	/*******释放内存*********/
	free(h_a);
	free(h_b);
	free(h_c);
	free(h_bt);
	cudaFree(d_a);
	cudaFree(d_b);
	cudaFree(d_c);
	cudaFree(d_bt);
	cudaFree(d_data);

	std::cout << "运行结束" << std::endl;
	return 0;
}