CUDA Matrix Transpose using Shared Memory 使用cuda共享内存优化矩阵转置

常常我们使用shared memory，为了解决以下两种问题：
1）同一个全局内存既被读取，又被写入。为了避免写入覆盖，于是放shared memory里暂存，所有数据处理好后，最后一起写入全局内存。
2）避免重复读取全局内存影响性能，所以一开始把数据存到shared memory，以便之后多次读取

本文将解决第3种问题：
3）为了全局内存的合并读取和写入，提高性能。
见如下矩阵转置的例子：

第一步，如图1，本thread将矩阵中对应位置的数据写到shared memory里面对应的位置
第二步，如图2.1，本thread将shared memory中对应的数据，转移到输出矩阵的转置的位置
如果按上面方法，使用shared memory并没有解决问题: 第一步，一个warp的threads对全局内存一行读取，进行的是合并读取，没毛病。第二步，对全局内存的一列写入，是离散的读取，效率很低，这个shared memory不仅没有解决问题，还增加了额外的copy时间。

所以我们要改变第二步，见图2.2：
为了读输出矩阵的写入是行写入，我们需要在第二步的时候对shared memory的数据进行重新index。
比如，在3列4行的block(我们申请的shared memory也是3列4行)里, 第7个thread的坐标是(2,1)，那么也找到shared memory的第7个位置，也就是(2,1)，所以它在第一步将输入数据copy到shared memory的(2,1)位置；第二步，重新index,方法是，转置之后的shared memory 找第7个位置，也就是(2,0)，所以第二步是将shared memory里(2,0)位置的数据copy到对应的输出矩阵中

代码如下：

#include <iostream>

using namespace std;
//    N
// M
#define N 8
#define M 4
#define BLOCKSIZE_X 2
#define BLOCKSIZE_Y 2

__global__ void launch_kernel(int *matrix, int *out){
	__shared__ int shmem[BLOCKSIZE_Y][BLOCKSIZE_X];

	// x,  y is the global thread location
	const unsigned  int x = blockIdx.x * blockDim.x + threadIdx.x;
	const unsigned int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(x>=N || y>=M) return;

	// threadIdx.y , threadIdx.x is the local shared memory index
	shmem[threadIdx.y][threadIdx.x] = matrix[y * N + x];
	__syncthreads();


	// find shmem[row][col] mapped to current thread
	const unsigned int tid = blockDim.x * threadIdx.y + threadIdx.x;
	const unsigned int col = tid / blockDim.y;
	const unsigned int row = tid % blockDim.y;

	// find the global index of (sx, sy) mapped to the local shmem[row][col]
	const unsigned int sx = blockIdx.x * blockDim.x + col;
	const unsigned int sy = blockIdx.y * blockDim.y + row;

	// then copy shmem[row][col] to out matrix (sy, sx)
	out[sx*M + sy] = shmem[row][col];

}

int main(){
	int matrix_h_in[M][N], *matrix_d_in;
	cout << "Original:" <<endl;
	for(int i=0; i<M; i++){
		for(int j=0; j<N; j++){
			matrix_h_in[i][j] =  rand()%10;
			cout << matrix_h_in[i][j] << "\t";
		}
		cout << endl;
	}
	cudaMalloc(&matrix_d_in, M*N*sizeof(int));
	cudaMemcpy(matrix_d_in, matrix_h_in, M*N*sizeof(int), cudaMemcpyHostToDevice);

	int matrix_h_out[N][M], *matrix_d_out;
	cudaMalloc(&matrix_d_out, N*M*sizeof(int));
	dim3 block(BLOCKSIZE_X,BLOCKSIZE_Y,1);
	dim3 grid((N + block.x-1)/block.x, (M + block.y-1)/block.y, 1);
	launch_kernel<<<grid,block>>>(matrix_d_in, matrix_d_out);
	cudaDeviceSynchronize();
	cudaMemcpy(matrix_h_out, matrix_d_out, N*M*sizeof(int), cudaMemcpyDeviceToHost);

	cout << "Transposed:" <<endl;
        for(int i=0; i<N; i++){
		for(int j=0; j<M; j++){
			cout << matrix_h_out[i][j] << "\t";
		}
		cout << endl;
	}

	cudaFree(matrix_d_in);
	cudaFree(matrix_d_out);
	

}

效果：

$ nvcc cudaMatrixTransposeSharedMemory.cu 
$ ./a.out 
Original:
3	6	7	5	3	5	6	2	
9	1	2	7	0	9	3	6	
0	6	2	6	1	8	7	9	
2	0	2	3	7	5	9	2	
Transposed:
3	9	0	2	
6	1	6	0	
7	2	2	2	
5	7	6	3	
3	0	1	7	
5	9	8	5	
6	3	7	9	
2	6	9	2