CUDA:矩阵转置的GPU实现(Share Memory)_cuda 矩阵转置

本文参加2022CUDA on Platform线上训练营学习笔记

欢迎各位大犇提意见

一、矩阵转置(Matrix Transpose)基础

上图中将m * n的矩阵A通过矩阵转置变成了n * m的 A^T,简单来讲矩阵转置即为将原始矩阵的第一行转置为目标矩阵的第一列，以此类推,相信基础扎实的你简单地看看CPU端的代码就能理解

二、矩阵转置的CPU端实现

__host__ void cpu_transpose(int *matrix,int *tr_matrix,int m,int n) {
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < m; j++) {
			tr_matrix[i * m + j] = matrix[j * n + i];
		}
	}
	return;
}
1
2
3
4
5
6
7
8

定义一个名为cpu_transpose的函数，将矩阵matrix转置为矩阵tr_matrix，通过观察代码不难发现tr_matrix[i][j]=matrix[j][i],这里需要注意到的是坐标的转换,转置后的矩阵行数和列数发生变换，留意m和n不要乘错了。

原始矩阵:

CPU端执行结果:

三、矩阵转置的GPU端实现(share Memory)

1、核函数的编写

GPU端的实现与CPU端类似，首先根据各个线程的index(索引)计算出当前线程在原始矩阵中的位置row和col,在原始矩阵中的row行，col列

	int row = blockDim.y * blockIdx.y + threadIdx.y;
	int col = blockDim.x * blockIdx.x + threadIdx.x;
1
2

下边我们申请同一个block中的线程可以访问的shared Memory

__shared__ int smem_matrix[BLOCK_SIZE][BLOCK_SIZE];
1

在GPU中申请了一块名为smem_matrix大小为sizeof(int)*BLOCK_SIZE^2的共享内存，在执行赋值操作之前将当前block中的线程需要访问到的数据从Global_Memory中复制到share_Memory

smem_matrix[threadIdx.y][threadIdx.x] = row < m&& col < n ? matrix[row*n+col] : 0;
1

赋值时需要注意的是:由于我们为内核函数设置执行配置的时候通常会向上取整，会申请多于实际需求的线程数，所以在我们赋值之前需要判断当前线程的坐标是否是需求坐标，以此来防止访问matrix时row*n+col成为野指针，对我们的数据造成重大的危害
有了同一个block中的线程申请一个share Memory的概念后，需要做的是同步同一个BLock中的线程

__syncthreads();
1

通过上边一系列的操作，我们就可以开始真正的转置操作了,需要注意的是，我们已经把线程所需的数据赋值到share Memory当中,所以我们在赋值时只需调用smem_matrix,同样,赋值操作之前，我们需要判断当前的坐标是否实际有效

	if(blockIdx.x * blockDim.x + threadIdx.y < n && threadIdx.x + blockIdx.y * blockDim.x < m)
	tr_matrix[threadIdx.x+blockIdx.y*blockDim.x+m*(blockIdx.x*blockDim.x+threadIdx.y)] = smem_matrix[threadIdx.x][threadIdx.y];
1
2

上述分析使我们获得了完整的GPU代码

__global__ void cuda_transpose(int *matrix,int *tr_matrix,int m,int n) {
	int row = blockDim.y * blockIdx.y + threadIdx.y;
	int col = blockDim.x * blockIdx.x + threadIdx.x;
	__shared__ int smem_matrix[BLOCK_SIZE][BLOCK_SIZE];
	smem_matrix[threadIdx.y][threadIdx.x] = row < m&& col < n ? matrix[row*n+col] : 0;
	__syncthreads();
	if(blockIdx.x * blockDim.x + threadIdx.y < n && threadIdx.x + blockIdx.y * blockDim.x < m)
	tr_matrix[threadIdx.x+blockIdx.y*blockDim.x+m*(blockIdx.x*blockDim.x+threadIdx.y)] = smem_matrix[threadIdx.x][threadIdx.y];
	return;
}

1
2
3
4
5
6
7
8
9
10
11

2、核函数的启动

在设备端申请两个指针并为其分配内存

	int* d_matrix, *dtr_matrix;
	cudaMalloc((void**)&d_matrix, sizeof(int) * m * n);
	cudaMalloc((void**)&dtr_matrix, sizeof(int) * m * n);
1
2
3

手动将matrix中的数据通过Pcie复制到设备端的Global Memory当中

cudaMemcpy(d_matrix, matrix, sizeof(int) * m * n, cudaMemcpyHostToDevice);
1

核函数执行设置的设定,一个warp通常为32个线程所以我们一个Block中的线程数最好设置为32的整数倍，从此提高使用率，有效防止inactive code的出现

dim3 block = { BLOCK_SIZE,BLOCK_SIZE,1 }; //BLOCK_SIZE = 16
1

gridDim的设置最需关注的就是申请的线程能够有效的覆盖真个矩阵，宁可多申请，通过核函数中的if屏蔽，也不少申请，导致计算的缺失，所以我们在计算中采用向上取整的方法
需要注意的使 dim3 类型中的三个成员都是要求unsigned int 类型的所以我们在前面添加(unsigned int)来强制将我们的计算结果转换为无符号

dim3 gird = { (unsigned int)(n - 1 + BLOCK_SIZE) / BLOCK_SIZE, (unsigned int)(m - 1 + BLOCK_SIZE) / BLOCK_SIZE,1 };
1

核函数启动！

cuda_transpose << < gird , block  >> > (d_matrix, dtr_matrix, m, n);
1

3、核函数性能计数

在CUDA中有一种特殊的类型cudaEvent_t，可以帮助我们记录核函数的执行信息

	cudaEvent_t kernel_start;
	cudaEvent_t kernel_end;

	cudaEventCreate(&kernel_start);
	cudaEventCreate(&kernel_end);
1
2
3
4
5

kernel_start用于记录核函数开始执行时的信息，kernel_end用来记录核函数运行结束时的信息，这里使用到了两个函数cudaEventQuery(kernel_start);,cudaEventSynchronize(kernel_end);,前者是非阻塞的，只要执行到就直接记录，后者是阻塞式的，需要前面的执行完毕才能运行,具体的性能计数函数如下

通过简单的逻辑组合，就可以得到核函数的实际运行时间，具体代码如下

	cudaEventCreate(&kernel_start);
	cudaEventCreate(&kernel_end);
	cudaEventRecord(kernel_start);
	cudaEventQuery(kernel_start);
	cuda_transpose << < gird , block  >> > (d_matrix, dtr_matrix, m, n);
	cudaEventRecord(kernel_end);
	cudaEventSynchronize(kernel_end);
	float ms;
	cudaEventElapsedTime(&ms, kernel_start, kernel_end);
1
2
3
4
5
6
7
8
9

四、代码参考

#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <iostream>
#define BLOCK_SIZE 32
using namespace std;


__global__ void cuda_transpose(int *matrix,int *tr_matrix,int m,int n) {
	int row = blockDim.y * blockIdx.y + threadIdx.y;
	int col = blockDim.x * blockIdx.x + threadIdx.x;
	__shared__ int smem_matrix[BLOCK_SIZE][BLOCK_SIZE];
	smem_matrix[threadIdx.y][threadIdx.x] = row < m&& col < n ? matrix[row*n+col] : 0;
	__syncthreads();
	if(blockIdx.x * blockDim.x + threadIdx.y < n && threadIdx.x + blockIdx.y * blockDim.x < m)
	tr_matrix[threadIdx.x+blockIdx.y*blockDim.x+m*(blockIdx.x*blockDim.x+threadIdx.y)] = smem_matrix[threadIdx.x][threadIdx.y];
	return;
}

__host__ void cpu_transpose(int *matrix,int *tr_matrix,int m,int n) {
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < m; j++) {
			tr_matrix[i * m + j] = matrix[j * n + i];
		}
	}
	return;
}

__host__ void init_matrix(int* matrix,int m,int n) {
	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) {
			matrix[i*n+j] = rand();
		}
	}
}

void print(int*, string,int,int);
bool check(int*, int*, int, int);

int main() {
	int m = 1111;
	int n = 113;
	int *matrix;
	cudaMallocHost((void**)&matrix, sizeof(int) * m * n);
	init_matrix(matrix,m,n);
	//print(matrix, "init matrix", m, n);


	int* htr_matrix;
	cudaMallocHost((void**)&htr_matrix, sizeof(int) * m * n);
	cpu_transpose(matrix, htr_matrix, m, n);
	//print(htr_matrix, "CPU", n, m);
	//将CPU端执行的结果存放在htr_matrix中 

	int* d_matrix, *dtr_matrix;
	cudaMalloc((void**)&d_matrix, sizeof(int) * m * n);
	cudaMalloc((void**)&dtr_matrix, sizeof(int) * m * n);
	cudaMemcpy(d_matrix, matrix, sizeof(int) * m * n, cudaMemcpyHostToDevice);
	dim3 gird = { (unsigned int)(n - 1 + BLOCK_SIZE) / BLOCK_SIZE, (unsigned int)(m - 1 + BLOCK_SIZE) / BLOCK_SIZE,1 };
	dim3 block = { BLOCK_SIZE,BLOCK_SIZE,1 };


	cudaEvent_t kernel_start;
	cudaEvent_t kernel_end;

	cudaEventCreate(&kernel_start);
	cudaEventCreate(&kernel_end);
	cudaEventRecord(kernel_start);
	cudaEventQuery(kernel_start);
	cuda_transpose << < gird , block  >> > (d_matrix, dtr_matrix, m, n);
	cudaEventRecord(kernel_end);
	cudaEventSynchronize(kernel_end);
	float ms;
	cudaEventElapsedTime(&ms, kernel_start, kernel_end);


	int* hdtr_matrix;
	cudaMallocHost((void**)&hdtr_matrix, sizeof(int) * m * n);
	cudaMemcpy(hdtr_matrix, dtr_matrix, sizeof(int) * m * n, cudaMemcpyDeviceToDevice);
	//print(hdtr_matrix, "GPU", n, m);
	
	if (check(hdtr_matrix, htr_matrix, n, m)) {
		cout << "pass\n";
	}
	else {
		cout << "error\n";
	}
	
	printf("GPU time is : %f \n", ms);

	cudaFree(hdtr_matrix);
	cudaFree(dtr_matrix);
	cudaFree(matrix);
	cudaFree(htr_matrix);
	cudaFree(d_matrix);
	return 0;
}


void print(int* a, string name,int m,int n) {
	cout << "NAME : " << name << endl;
	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) {
			printf("%6d ", a[i * n + j]);
		}
		printf("\n");
	}
}

bool check(int* a, int* b, int m, int n) {
	bool check_flag = true;
	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) {
			if (a[i * n + j] != b[i * n + j]) {
				return false;
			}
		}
	}
	return check_flag;
}


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124

执行结果如图

五、实践心得

本次实践通过GPU端中的share Memory对核函数运行时的读写问题做了优化，当线程与线程之间为连续读写时，global Memory的效率是比较高的，不使用share Memory时，使用GPU进行矩阵转置会出现两难问题（1.读row-major 写col-major，2写col-major 读row-major），而在share Memory中row-major和col-major的效率几乎相同，很好地解决了global memory上的问题，在编写过程中，需要注意的是，要顺着global memory写，首先保证global memory读写时是row-major，以达到最高的优化效率。
遇到的最大问题是，边界的判断问题，GPU转置过程中，由于要保证global memory 是 row-major，所以坐标不像是CPU端中的简单调换，具体表现为（在对share 数字赋值时该线程无意义，而在写global操作中该线程有意义），所以在__syncthreads();后需要判断当前线程是否有意义

鄙人第一次写实操博客，有建议必洗耳恭听
再次感谢伟大的NV 开发者社区