gpu cuda 矩阵乘法

欢迎关注更多精彩
关注我，学习常用算法与数据结构，一题多解，降维打击。

问题描述

给定2个二数组表示矩阵，利用gpu相乘并返回结果。

cpu 算法

void cpu_matrix_multi() {
	for (int y = 0; y < M; ++y) {
		for (int x = 0; x < M; ++x) {
			cpu_result[y][x] = 0;
			for (int k = 0; k < N; ++k) {
				cpu_result[y][x] += matrixa[y][k] * matrixb[k][x];
			}
		}
	}
}

cpu耗时

耗时786毫秒

gpu 算法-直接法

代码

__global__ void gpu_matrix_multi(int ina[M][N], int inb[N][M], int out[M][M]) {
	int x = threadIdx.x + blockDim.x*blockIdx.x;
	int y = threadIdx.y + blockDim.y*blockIdx.y;

	if (x < M && y < M) {
		out[y][x] = 0;
		for (int k = 0; k < N; ++k) {
			out[y][x] += ina[y][k] * inb[k][x];
		}
	}
}

耗时分析

耗时6.504毫秒，提升了100多倍。

算法缺点

在gpu中数组是放在全局内存里，全局内存访问比较慢（相较于共享内存）。如果是连续访问则可以进行合并访存，效率上可以作一些弥补，如果一直随机访问效率会打折扣。

分析上述代码。

对于初始矩阵A[M][K], B[K][N]相乘，得到C[M][N].

那么对于A[i][j] 单个素要被访问N次，对于B[i][j]单个元素要被访问M次，每次都要从全局内存中访问，非常慢。

对于C[i][j] 也要被写入K次。

gpu 算法分块法（块内顺序访问）

算法思路

共享内存访问效率比全局内存快速1个数量级。

但是共享内存大小有限制。

gpu中1个block内存所有线程可以共同访问一块共享内存。

block最多有1024个线程。

我们可以把矩阵分成多个（n*m）块，每个块32*32大小。

每个block可以分批进行乘法聚合。

聚合1个块具体过程如下：

步骤1： 以计算out[6]块为例，先将in1, in2分块，取到out[6]对应的相关数据块，分成in1[1,2,3], int2[1,2,3]

步骤2： 把in1[1,2,3], in2[1,2,3] 依次加载到共享内存，相乘后加到out[6]块。

即out[6] = in1[1]*in2[1] + in1[2]*in2[2] + in1[3]*in2[3]

这样在访问in, out时都做到了访存复用，效率得到提升。

代码

#include "cuda_runtime.h"
#include <math.h>
#include<vector>
#include "device_launch_parameters.h"

#include <stdio.h>

using namespace std;

#define BLOCK_SIZE 32
#define M 1000
#define N 500

__managed__ int matrixa[M][N];
__managed__ int matrixb[N][M];
__managed__ int gpu_result[M][M];
__managed__ int cpu_result[M][M];

__global__ void gpu_matrix_multi(int ina[M][N], int inb[N][M], int out[M][M]) {
	int x = threadIdx.x + blockDim.x*blockIdx.x;
	int y = threadIdx.y + blockDim.y*blockIdx.y;

	if (x < M && y < M) {
		out[y][x] = 0;
		for (int k = 0; k < N; ++k) {
			out[y][x] += ina[y][k] * inb[k][x];
		}
	}
}

__global__ void gpu_shared_matrix_multi(int ina[M][N], int inb[N][M], int out[M][M]) {
	int x = threadIdx.x + blockDim.x*blockIdx.x;
	int y = threadIdx.y + blockDim.y*blockIdx.y;

	__shared__ int sa[BLOCK_SIZE + 1][BLOCK_SIZE + 1];
	__shared__ int sb[BLOCK_SIZE + 1][BLOCK_SIZE + 1];
	int temp = 0;
	for (int step = 0; step <= N / BLOCK_SIZE; step++) {
		int step_x = step * BLOCK_SIZE + threadIdx.x;
		if (y < M && step_x < N)
			sa[threadIdx.y][threadIdx.x] = ina[y][step_x];
		else
			sa[threadIdx.y][threadIdx.x] = 0;

		int step_y = step * BLOCK_SIZE + threadIdx.y;

		if (step_y < M && x < N) 
			sb[threadIdx.y][threadIdx.x] = inb[step_y][x];
		else
			sb[threadIdx.y][threadIdx.x] = 0;

		__syncthreads();

		for (int k = 0; k < BLOCK_SIZE; ++k) {
			temp += sa[threadIdx.y][k] * sb[k][threadIdx.x];
		}

		__syncthreads();
	}
	if(y<M&&x<M)out[y][x] = temp;
}

void cpu_matrix_multi() {
	for (int y = 0; y < M; ++y) {
		for (int x = 0; x < M; ++x) {
			cpu_result[y][x] = 0;
			for (int k = 0; k < N; ++k) {
				cpu_result[y][x] += matrixa[y][k] * matrixb[k][x];
			}
		}
	}
}


int main()
{
	for (int y = 0; y < N; ++y) {
		for (int x = 0; x < M; ++x) {
			matrixa[y][x] = x + y * M;
			matrixb[x][y] = x + y * M;
		}
	}

	cudaEvent_t start, stop_gpu, stop_cpu;
	cudaEventCreate(&start);
	cudaEventCreate(&stop_gpu);
	cudaEventCreate(&stop_cpu);

	cudaEventRecord(start);
	cudaEventSynchronize(start);

	dim3 dimGrid((M+BLOCK_SIZE-1)/BLOCK_SIZE,(M+BLOCK_SIZE-1)/BLOCK_SIZE);
	dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);

		gpu_shared_matrix_multi<<<dimGrid,dimBlock>>>(matrixa, matrixb, gpu_result);
		//gpu_matrix_multi<<<dimGrid,dimBlock>>>(matrixa, matrixb, gpu_result);
		cudaDeviceSynchronize();

	cudaEventRecord(stop_gpu);
	cudaEventSynchronize(stop_gpu);

	cpu_matrix_multi();
	cudaEventRecord(stop_cpu);
	cudaEventSynchronize(stop_cpu);

	float time_cpu, time_gpu;
	cudaEventElapsedTime(&time_gpu, start, stop_gpu);
	cudaEventElapsedTime(&time_cpu, stop_gpu, stop_cpu);

	bool errors = false;
	for (int y = 0; y < M; ++y) {
		for (int x = 0; x < M; ++x) {
			if (cpu_result[y][x] != gpu_result[y][x]) errors = true;
		}
	}

	printf("result: %s\n", errors? "fault":"pass");
	printf("CPU time: %.3f\nGPU time: %.3f\n", time_cpu, time_gpu);

	cudaEventDestroy(start);
	cudaEventDestroy(stop_gpu);
	cudaEventDestroy(stop_cpu);

    return 0;
}

耗时分析

耗时4.702毫秒，比直接法提升了20%多。

---

本人码农，希望通过自己的分享，让大家更容易学懂计算机知识。创作不易，帮忙点击公众号的链接。