本文通过编写CUDA核函数，以矩阵转置等5个问题为例，利用GPU众核对大规模问题的求解加速。
矩阵转置
算法流程
- 将需要转置的矩阵存储到GPU内存中。
- 在GPU上分配空间存储转置后的矩阵。
- 定义CUDA核函数来实现矩阵转置。该核函数应该使用线程块和线程格的概念来处理矩阵中的所有元素。在每个线程块中，线程可以使用共享内存来处理数据。最后，利用全局内存将结果写回到GPU。
- 调用CUDA核函数以执行矩阵转置。
- 将转置后的矩阵从GPU内存复制到主机内存中。
- 释放GPU内存
代码实现
在传统代码的基础上，利用共享内存优化访问全局内存的方式，同时使用padding操作，避免bank冲突。
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
|
const int N = 10000;
const int TILE_DIM = 32;
const int grid_size_x = (N + TILE_DIM - 1) / TILE_DIM;
const int grid_size_y = grid_size_x;
const dim3 block_size(TILE_DIM, TILE_DIM);
const dim3 grid_size(grid_size_x, grid_size_y);
__global__ void transpose(const float *A, float *B, const int N){
__shared__ float S[TILE_DIM][TILE_DIM+1];
int bx = blockIdx.x * TILE_DIM;
int by = blockIdx.y * TILE_DIM;
//顺序读
int nx1 = bx + threadIdx.x;
int ny1 = by + threadIdx.y;
if (nx1 < N && ny1 < N)
S[threadIdx.y][threadIdx.x] = A[ny1 * N + nx1];
__syncthreads();
//顺序写
int nx2 = bx + threadIdx.y;
int ny2 = by + threadIdx.x;
if (nx2 < N && ny2 < N)
B[nx2 * N + ny2] = S[threadIdx.x][threadIdx.y];
}
int main(){
int nn = N * N;
int nBytes = nn * sizeof(float);
float a = 1.1;
float b = 2.2;
//cpu空间分配
float* A = (float*)malloc(nBytes);
float* B = (float*)malloc(nBytes);
//矩阵初始化
for (int i = 1; i <= nn; i++) {
if (i % 2 == 0)
A[i-1] = a;
else
A[i-1] = b;
}
//gpu空间分配（需要先定义）
float* d_A, * d_B;
cudaMalloc((void**)&d_A, nBytes);
cudaMalloc((void**)&d_B, nBytes);
cudaMemcpy(d_A, A, nBytes, cudaMemcpyHostToDevice);
transpose<<<grid_size, block_size >>> (d_A, d_B, N);
cudaMemcpy(B, d_B, nBytes, cudaMemcpyDeviceToHost);
free(A);
free(B);
cudaFree(d_A);
cudaFree(d_B);
return 0;
}
|
矩阵相乘
核函数设计
共享内存的使用：应该将两个原始矩阵的一部分抽取出来，当作瓦片，放入共享内存，减少对全局内存的访问，从而提高程序性能。
非方阵处理：非方阵的两个矩阵相乘，需要调整共享内存中瓦片的尺寸，同时需要添加判断语句，防止发生越界错误。
代码实现
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
|
__global__ void matrixmul(const float *A, const float *B, float *C,const int N,const int M,const int K){ //N*M x M*K
__shared__ float A_S[TILE_DIM][TILE_DIM];
__shared__ float B_S[TILE_DIM][TILE_DIM];
int tx = threadIdx.x;
int ty = threadIdx.y;
int bx = blockIdx.x;
int by = blockIdx.y;
int row = by * TILE_DIM + ty;
int col = bx * TILE_DIM + tx;
float value = 0;
for (int ph = 0; ph < M / TILE_DIM+1; ph++) {
if (row < N && ph * TILE_DIM + tx < M)
A_S[ty][tx] = A[row * M + ph * TILE_DIM + tx];
else
A_S[ty][tx] = 0.0;
if (col < K && ph * TILE_DIM + ty < M)
B_S[ty][tx] = B[(ph * TILE_DIM + ty) * K + col];
else
B_S[ty][tx] = 0.0;
__syncthreads();
for (int k = 0; k < TILE_DIM; k++)
value += A_S[ty][k] * B_S[k][tx];
__syncthreads();
}
if (row < N && col < K)
C[row*K+col]=value;
}
|
统计直方图
共享内存的使用：在共享内存开辟临时数组，减少对全局内存的访问，从而提高程序性能。
跨网格循环：数据量很大，一个线程要处理多个数据。
原子操作：无论是共享内存还是全局内存的结果，都需要使用原子操作，避免冲突
1
2
3
4
5
6
7
8
9
10
11
12
13
|
__global__ void histo_kernel(unsigned char *buffer, long size, unsigned int *histo){
__shared__ unsigned int temp[256];
temp[threadIdx.x] = 0;
__syncthreads();
int i = threadIdx.x + blockIdx.x * blockDim.x;
int offset = blockDim.x *gridDim.x;
while (i<size){
atomicAdd(&temp[buffer[i]], 1);
i += offset;
}
__syncthreads();
atomicAdd(&(histo[threadIdx.x]), temp[threadIdx.x]);
}
|
规约求和
跨网格循环：对于海量数据(>100w)，即使是GPU众核，也做不到这么多的线程，因此一个线程要处理多个数据。
共享内存优化：将每个线程对多个元素求和的结果存于共享内存，可以减少对全局内存的访问。
交错配对法：利用交错配对法进行规约，可以缓解一部分的线程束分化现象。
原子操作:将结果累加到全局内存中，可能会发生冲突，需要使用原子操作。
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
|
__global__ void _sum_gpu(int *input, int count, int *output)
{
__shared__ int sum_per_block[BLOCK_SIZE];
int temp = 0;
for (int idx = threadIdx.x + blockDim.x * blockIdx.x;
idx < count; idx += gridDim.x * blockDim.x
)
temp += input[idx];
sum_per_block[threadIdx.x] = temp;
__syncthreads();
//**********shared memory summation stage***********
for (int length = BLOCK_SIZE / 2; length >= 1; length /= 2)
{
int double_kill = -1;
if (threadIdx.x < length)
double_kill = sum_per_block[threadIdx.x] + sum_per_block[threadIdx.x + length];
__syncthreads(); //why we need two __syncthreads() here, and,
if (threadIdx.x < length)
sum_per_block[threadIdx.x] = double_kill;
__syncthreads();
} //the per-block partial sum is sum_per_block[0]
if (threadIdx.x == 0) atomicAdd(output, sum_per_block[0]);
|
TOP K 问题
TOP K问题是指在一组数据中，找到前K个最大或最小的元素。利用CUDA规约计算可以高效地解决TOP K问题。
以下是利用CUDA规约计算来实现排序和选择前K个最大/最小元素的详细步骤：
- 定义一个二元组类型，包含值和索引，用于存储原始数据及其索引（为了在排序后恢复原始数据）
- 对原始数据进行遍历，将数据存储到二元组中
- 对二元组进行归约操作，得到前K个最大/最小值的索引
- 恢复原始数据并按照索引排序，得到前K个最大/最小值
具体实现流程如下：
- 将数据复制到GPU显存中
float *d_data; cudaMemcpy(d_data, h_data, size, cudaMemcpyHostToDevice);
- 将数据存储到二元组中
1
2
3
4
5
|
typedef struct { float value; int index; } Tuple;
Tuple *d_tuples;
int threadsPerBlock = 256;
int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;
initializeTuples<<<blocksPerGrid, threadsPerBlock>>>(d_data, d_tuples, n);
|
- 对二元组进行归约操作，得到前K个最大/最小值的索引
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
|
int *d_indices;
kReduceKernel<<<blocksPerGrid, threadsPerBlock>>>(d_tuples, d_indices, n, k);
__global__ void kReduceKernel(Tuple *input, int *output, int n, int k) {
extern __shared__ Tuple shared[];
int tid = threadIdx.x;
int i = blockIdx.x * blockDim.x + threadIdx.x;
shared[tid] = (i < n) ? input[i] : Tuple{0, 0};
__syncthreads();
for (int s = blockDim.x / 2; s > 0; s >>= 1) {
if (tid < s)
shared[tid] = (shared[tid].value > shared[tid + s].value) ? shared[tid] : shared[tid + s];
__syncthreads();
}
if (tid == 0)
output[blockIdx.x] = shared[0].index;
}
|
- 在CPU中恢复原始数据并按照索引排序，得到前K个最大/最小值
1
2
3
4
5
|
cudaMemcpy(h_indices, d_indices, size, cudaMemcpyDeviceToHost);
for (int i = 0; i < k; ++i) {
int index = h_indices[i];
h_result[i] = h_data[index]; }
std::sort(h_result, h_result + k);
|
Tim
支付宝
贝宝
微信