并行编程实验二
实验目的
- 掌握CUDA网格和线程块的设计
- 实践CUDA运行时API的使用
- 通过编写核函数，掌握利用GPU众核对大规模问题的求解加速的方法
- 体会CPU多核与GPU众核这两种不同体系结构带来的程序设计上的差异
实验内容
矩阵转置
算法流程
- 将需要转置的矩阵存储到GPU内存中。
- 在GPU上分配空间存储转置后的矩阵。
- 定义CUDA核函数来实现矩阵转置。该核函数应该使用线程块和线程格的概念来处理矩阵中的所有元素。在每个线程块中，线程可以使用共享内存来处理数据。最后，利用全局内存将结果写回到GPU。
- 调用CUDA核函数以执行矩阵转置。
- 将转置后的矩阵从GPU内存复制到主机内存中。
- 释放GPU内存
代码实现
在传统代码的基础上，利用共享内存优化访问全局内存的方式，同时使用padding操作，避免bank冲突。
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;
}矩阵相乘
核函数设计
共享内存的使用：应该将两个原始矩阵的一部分抽取出来，当作瓦片，放入共享内存，减少对全局内存的访问，从而提高程序性能。
非方阵处理：非方阵的两个矩阵相乘，需要调整共享内存中瓦片的尺寸，同时需要添加判断语句，防止发生越界错误。
代码实现
__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;
}统计直方图
共享内存的使用：在共享内存开辟临时数组，减少对全局内存的访问，从而提高程序性能。
跨网格循环：数据量很大，一个线程要处理多个数据。
原子操作：无论是共享内存还是全局内存的结果，都需要使用原子操作，避免冲突
__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众核，也做不到这么多的线程，因此一个线程要处理多个数据。
共享内存优化：将每个线程对多个元素求和的结果存于共享内存，可以减少对全局内存的访问。
交错配对法：利用交错配对法进行规约，可以缓解一部分的线程束分化现象。
原子操作:将结果累加到全局内存中，可能会发生冲突，需要使用原子操作。
__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); - 将数据存储到二元组中
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); ``` 3. 对二元组进行归约操作，得到前K个最大/最小值的索引 ```cpp 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个最大/最小值
完成以上步骤后，就可以得到排序后的前K个最大/最小值了。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);