高性能计算面试笔记

CUDA kernel

elementwise

逐元素操作,难度不高,考频低,本文以add为例。

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
//	Naive GPU
template <typename ValueType, typename IndexType>
__global__ void addKernel_naive(ValueType *A, ValueType *B, ValueType *C, IndexType N) {
IndexType idx = blockDim.x * blockIdx.x + threadIdx.x; // global index
if(idx < N) C[idx] = A[idx] + B[idx];
}

// 以 ValueType 为 float 为例进行向量化访存
__global__ void addKernel_float4(float *A, float *B, float *C, int N) {
int idx = (blockDim.x * blockIdx.x + threadIdx.x) * 4; // 4 floats a time idx * 4
if (idx < N) {
float4 tmpa = reinterpret_cast<float4*>(A)[idx / 4];
float4 tmpb = reinterpret_cast<float4*>(B)[idx / 4];
float4 tmpc;

tmpc.x = tmpa.x + tmpb.x;
tmpc.y = tmpa.y + tmpb.y;
tmpc.z = tmpa.z + tmpb.z;
tmpc.w = tmpa.w + tmpb.w;

reinterpret_cast<float4*>(C)[idx / 4] = tmpc;
}
}

// 向量化访存让一个线程处理了更多数据,我们可以进一步加大stride,让一个线程处理更多个float4
__global__ void addKernel_float4_stride(float *A, float *B, float *C, int N) {
int start_idx = blockDim.x * blockIdx.x + threadIdx.x;
int stride = blockDim.x * gridDim.x;

for (int i = start_idx; i * 4 < N; i += stride) {
int base_idx = i * 4;

// 如果剩下的元素足够 4 个
if (base_idx + 3 < N) {
float4 tmp_a = reinterpret_cast<float4*>(A)[i];
float4 tmp_b = reinterpret_cast<float4*>(B)[i];
float4 res;
res.x = tmp_a.x + tmp_b.x;
res.y = tmp_a.y + tmp_b.y;
res.z = tmp_a.z + tmp_b.z;
res.w = tmp_a.w + tmp_b.w;
reinterpret_cast<float4*>(C)[i] = res;
} else {
// 最后的尾部处理:逐个处理剩下的 1~3 个元素
for (int j = base_idx; j < N; ++j) {
C[j] = A[j] + B[j];
}
}
}
}

reduceSum

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
__global__ void reduceSum(float *d_x, float *d_y, int N) {
__shared__ float sy[32]; // 一个warp的大小

int idx = blockDim.x * blockIdx.x + threadIdx.x;
int warpIdx = threadIdx.x / warpSize; // 第几个warp
int laneIdx = threadIdx.x % warpSize; // 该warp中的第几个线程

float val = (idx < N) ? d_x[idx] : 0.0f;

for (int offset = warpSize >> 1; offset > 0; offset >>= 1) {
val += __shfl_down_sync(0xFFFFFFFF, val, offset);
}

if(laneIdx == 0) sy[warpIdx] = val; // 用每个warp的第一个线程将一个warp规约的结果存储在对应的smem位置
__syncthreads();

if(warpIdx == 0) {
int warpNum = blockDim.x / warpSize;
val = (laneIdx < warpNum) ? sy[laneIdx] : 0.0f;
for (int offset = warpSize >> 1; offset > 0; offset >>= 1) {
val += __shfl_down_sync(0xFFFFFFFF, val, offset);
}

if(laneIdx == 0) atomicAdd(d_y, val);
}
}

// float4打满带宽
__global__ void reduceSum_float4(float *d_x, float *d_y, int N) {
__shared__ float sy[32]; // 一个warp的大小

int idx = (blockDim.x * blockIdx.x + threadIdx.x) * 4;
int warpIdx = threadIdx.x / warpSize; // 第几个warp
int laneIdx = threadIdx.x % warpSize; // 该warp中的第几个线程
float val = 0.0f;

if (idx < N) {
float4 tmp_x = reinterpret_cast<float4*>(d_x)[idx / 4];
val += tmp_x.x;
val += tmp_x.y;
val += tmp_x.z;
val += tmp_x.w;
}

for (int offset = warpSize >> 1; offset > 0; offset >>= 1) {
val += __shfl_down_sync(0xFFFFFFFF, val, offset);
}

if(laneIdx == 0) sy[warpIdx] = val; // 用每个warp的第一个线程将一个warp规约的结果存储在对应的smem位置
__syncthreads();

if(warpIdx == 0) {
int warpNum = blockDim.x / warpSize;
val = (laneIdx < warpNum) ? sy[laneIdx] : 0.0f;
for (int offset = warpSize >> 1; offset > 0; offset >>= 1) {
val += __shfl_down_sync(0xFFFFFFFF, val, offset);
}

if(laneIdx == 0) atomicAdd(d_y, val);
}
}

// 一个线程处理更多数据
__global__ void reduceSum_float4_stride(float* d_x, float* d_y, const int N) {
__shared__ float s_y[32];

int stride = blockDim.x * gridDim.x * 4;
int idx = (blockDim.x * blockIdx.x + threadIdx.x) * 4;

int warpId = threadIdx.x / warpSize;
int laneId = threadIdx.x % warpSize;

float val = 0.0f;

for (int i = idx; i < N; i += stride) {
// 如果 N 不是 4 的倍数
if (i + 3 < N) {
float4 tmp_x = reinterpret_cast<float4*>(d_x)[i / 4];
val += tmp_x.x + tmp_x.y + tmp_x.z + tmp_x.w;
} else {
// 处理数组末尾不足 4 个元素的情况
for (int j = i; j < N; ++j) {
val += d_x[j];
}
}
}

for (int offset = warpSize >> 1; offset > 0; offset >>= 1) {
val += __shfl_down_sync(0xFFFFFFFF, val, offset);
}

if (laneId == 0) s_y[warpId] = val;
__syncthreads();

if (warpId == 0) {
int warpNum = blockDim.x / warpSize;
val = (laneId < warpNum) ? s_y[laneId] : 0.0f;

for (int offset = warpSize >> 1; offset > 0; offset >>= 1) {
val += __shfl_down_sync(0xFFFFFFFF, val, offset);
}

if (laneId == 0) atomicAdd(d_y, val);
}
}

softmax

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
// cpu version
void softmax(float *input, float *output, int size) {
float max = -INFINITY;
int i;

// 找到这组数据的最大值
for (i = 0; i < size; i++) {
if(input[i] > max) max = input[i];
}

float sum = 0.0f;
for (i = 0; i < size; i++) {
// 减去最大值防止溢出
output[i] = expf(input[i]-max);
sum += output[i];
}

for (i = 0; i < size; i++) {
output[i] /= sum;
}
}

// gpu version
__global__ void softmaxKernel(float *input, float *output, int cols, int rows) {
int row = blockDim.x * blockIdx.x + threadIdx.x;
if (row >= rows) return;
float max = -INFINITY;
int i;
for (i = 0; i < cols; i++) {
float tmp = input[row * cols + i];
if(tmp > max) max = tmp;
}
float sum = 0.f;
for (i = 0; i < cols; i++) {
float tmp = input[row * cols + i];
float tmp_exp = expf(tmp - max);
output[row * cols + i] = tmp_exp;
sum += tmp_exp;
}

for (i = 0; i < cols; i++) {
output[row * cols + i] /= sum;
}
}

sgemm

高频考察,需关注

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
//	cpu version
void Sgemm(int M, int N, int K, float alpha, float beta, float *A, float *B, float *C) {

for (int i = 0; i < M; i++) {
for (int j = 0; j < N; j++) {
float sum = 0.0f;
for (int k = 0; k < K; k++) {
sum += A[i * K + k] * B[k * N + j];
}
C[i * N + j] = alpha * sum + beta * C[i * N + j];
}
}
}

// gpu version
__global__ void SgemmKernel(int M, int N, int K, float alpha, float beta, float *A, float *B, float *C) {
int col = blockDim.x * blockIdx.x + threadIdx.x;
int row = blockDim.y * blockIdx.y + threadIdx.y;
if (row >= M || col >= N) return;

float sum = 0.f;
for (int k = 0; k < K; k++) {
sum += A[row * K + k] * B[k * N + col];
}
C[row * N + col] = alpha * sum + beta * C[row * N + col];
}

// share Memory
template<const int BLOCK_SIZE>
__global__ void SgemmKernel(int M, int N, int K, float alpha, float beta, float *A, float *B, float *C) {
int BM = BLOCK_SIZE;
int BN = BLOCK_SIZE;
int BK = BLOCK_SIZE;

int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x % BN;
int ty = threadIdx.y / BN;

// 偏移到当前线程的索引
A = &A[by * BM * K];
B = &B[bx * BN];
C = &C[by * BM * K + bx * BN];

__shared__ float As[BM * BK];
__shared__ float Bs[BK * BN];

float sum =0.f;

for (int i = 0; i < K; i+= BK) {
As[ty * BK + tx] = A[ty * K + tx];
Bs[ty * BN + tx] = B[ty * N + tx];
__syncthreads();

//偏移到下一步
A += BK;
B += BK * N;

for(int k = 0; k < BK; k++) {
sum += As[ty * BK + k] * Bs[k * BN + tx];
}
__syncthreads();

}
C[ty * N + tx] = alpha * sum + C[ty * N + tx];
}

transpose

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
//	0维和3维互换
void transpose(float *dst, float *src, int M, int N, int K, int W) {
// 原形状 M*N*K*W
int ostride0 = N * K * W;
int ostride1 = K * W;
int ostride2 = W;
// 转置后 W*N*K*M
int nstride0 = M * N * K;
int nstride1 = K * M;
int nstride2 = M;
for (int m = 0; m < M; m++) {
for (int n = 0; n < N; n++) {
for (int k = 0; k < K; k++) {
for (int w = 0; w < W; w++) {
int oldidx = m * ostride0 + n * ostride1 + k * ostride2 + w;
int newidx = w * nstride0 + n * nstride1 + k * nstride2 + m;
dst[newidx] = src[oldidx];
}
}
}
}
}

__global__ void transpose(float *dst, float *src, int M, int N, int K, int W) {
int m = blockDim.x * blockIdx.x + threadIdx.x;
int n = blockDim.y * blockIdx.y + threadIdx.y;
// 原形状 M*N*K*W
int ostride0 = N * K * W;
int ostride1 = K * W;
int ostride2 = W;
// 转置后 W*N*K*M
int nstride0 = M * N * K;
int nstride1 = K * M;
int nstride2 = M;

for (int k = 0; k < K; k++) {
for (int w = 0; w < W; w++) {
int oldidx = m * ostride0 + n * ostride1 + k * ostride2 + w;
int newidx = w * nstride0 + n * nstride1 + k * nstride2 + m;
dst[newidx] = src[oldidx];
}
}
}

conv2d

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
//	cpu version
void conv2d(const float *input, const float *kernel, float *output, int iH, int iW, int kH, int kW) {
int oH = iH - kH + 1;
int oW = iW - kW + 1;

for (int h = 0; h < oH; h++) {
for (int w = 0; w < oW; w++) {
float sum = 0.0f;
for (int ky = 0; ky < kH; ky++) {
for (int kx = 0; kx < kW; kx ++) {
sum += input[(h + ky) * iW + (w + kx)] * kernel[ky * kW + kx];
}
}
output[oW * h + w] = sum;
}
}
}

// gpu version
__global__ void conv2dKernel(const float *input, const float *kernel, float *output, int iH, int iW, int kH, int kW, int oH, int oW) {
int h = blockDim.y * blockIdx.y + threadIdx.y;
int w = blockDim.x * blockIdx.x + threadIdx.x;

if(h < oH && w < oW) {
float sum = 0.f;
for (int ky = 0; ky < kH; ky ++) {
for (int kx = 0; kx < kW; kx++) {
sum += input[(h + ky) * iW + (kx + w)] * kernel[ky * kW + kx];
}
}
output[oW * h + w] = sum;
}
}

sgemv

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
//cpu version
template <typename IndexType, typename ValueType>
void spmv_cpu_kernel(vector<IndexType> &row_offset,
vector<IndexType> &col_index,
vector<ValueType> &value,
vector<ValueType> &x,
vector<ValueType> &y,
IndexType row_num)
{
for(int i=0; i<row_num; i++){
ValueType res = 0;
IndexType num = row_offset[i+1] - row_offset[i];
for(int j=0; j<num; j++){
IndexType index = row_offset[i] + j;
res += value[index]*x[col_index[index]];
}
y[i] = res;
}
}

//gpu version
template <unsigned int WarpSize>
__device__ __forceinline__ float warpReduceSum(float sum) {
if (WarpSize >= 32)sum += __shfl_down_sync(0xffffffff, sum, 16); // 0-16, 1-17, 2-18, etc.
if (WarpSize >= 16)sum += __shfl_down_sync(0xffffffff, sum, 8);// 0-8, 1-9, 2-10, etc.
if (WarpSize >= 8)sum += __shfl_down_sync(0xffffffff, sum, 4);// 0-4, 1-5, 2-6, etc.
if (WarpSize >= 4)sum += __shfl_down_sync(0xffffffff, sum, 2);// 0-2, 1-3, 4-6, 5-7, etc.
if (WarpSize >= 2)sum += __shfl_down_sync(0xffffffff, sum, 1);// 0-1, 2-3, 4-5, etc.
return sum;
}

template <typename IndexType, typename ValueType, unsigned int VECTORS_PER_BLOCK, unsigned int THREADS_PER_VECTOR>
__global__ void My_spmv_csr_kernel(const IndexType row_num,
const IndexType * A_row_offset,
const IndexType * A_col_index,
const ValueType * A_value,
const ValueType * x,
ValueType * y)
{
const IndexType THREADS_PER_BLOCK = VECTORS_PER_BLOCK * THREADS_PER_VECTOR;
const IndexType thread_id = THREADS_PER_BLOCK * blockIdx.x + threadIdx.x; // global thread index
const IndexType thread_lane = threadIdx.x & (THREADS_PER_VECTOR - 1); // thread index within the vector
const IndexType row_id = thread_id / THREADS_PER_VECTOR; // global vector index

if(row_id < row_num){
const IndexType row_start = A_row_offset[row_id]; //same as: row_start = Ap[row];
const IndexType row_end = A_row_offset[row_id+1];

// initialize local sum
ValueType sum = 0;

// accumulate local sums
for(IndexType jj = row_start + thread_lane; jj < row_end; jj += THREADS_PER_VECTOR)
sum += A_value[jj] * x[ A_col_index[jj] ];

sum = warpReduceSum<THREADS_PER_VECTOR>(sum);
if (thread_lane == 0){
y[row_id] = sum;
}
}
}

Problems

高性能计算

1. SIMD 和 SIMT 的区别

simd 指“单指令多数据”,一般用于cpu,面向处理数据组,simt 则指“单指令多线程”,一般指面向gpu的线程组。

2. CUDA 线程组织

GRID(网格):一个核函数启动的所有线程的总和。

BLOCK(线程块):一个 grid 可以分成多个 block ,每个 block 有唯一标识可通过 blockIdx.x,blockIdx.y,blockIdx.z 访问,同一个 block 内的线程可以访问同一块共享内存(shared memory),可使用 syncthreads() 进行同步,block是资源分配的最小单位,gpu 以 block 为单位将线程分配到 sm ,一个 sm 可以放多个 block,一个 block 会固定到分配的 sm 到程序结束。

THREAD(线程):一个 block 里有多个 thread,thread 是最小的逻辑运行单位,可通过 threadIdx.x,threadIdx.y,threadIdx.z 进行访问。

3.warp的作用

GPU 内的线程被组织成每32个线程一组的 warp(线程束) ,因为在 SIMT 模式下,希望更多的线程能高效执行,如果每个线程都挨着启动统一指令,效率必然不高,如果所有线程同时启动同一指令,又丧失了灵活性,于是有了warp (线程束)。多个 warp 可以在 SM 上并行执行,同一个 warp 的线程执行同一条命令,同一个 warp 内的线程并行处理不同数据,从而在保证并行处理能力的同时提高了 gpu 的灵活性。

4.GPU SM

GPU 的线程调度主要由 SM (流式多处理器,核心计算单元,内部结构很多,接下来列举说明)中的 warp scheduler 负责。

Instruction cache(指令缓存):存放指令的缓存空间,减少 SM 从全局内存获取指令的延迟——不同的 warp 可以复用已经缓存好的指令,利用时间局部性(短时间内可能重复执行同一条指令)和空间局部性(程序的指令之间是连续的)。

Warp Scheduler(Warp调度器):决定每个时钟周期要执行的 warp ,当一个warp 由于内存访问或数据依赖而无法执行时,Warp Scheduler 会调度其他 warp 继续执行,确保 SM 内的计算单元尽可能不空闲,最大化资源利用率,同一个 warp 的线程同步执行同一条指令。

Dispatch Unit(指令分发单元):负责将 Warp Scheduler 选择的指令分发到合适的计算单元(如CUDA Cores、Tensor Cores)上的模块,Dispatch Unit 根据 Warp Scheduler 的决定,从 Instruction Cache 中获取已经调度好的指令,将指令发送到相应的硬件执行单元(如CUDA Core、ALU、Load/Store单元等)。它协调执行不同类型的指令,确保每个硬件单元处理适合其功能的指令。

SP(流处理器):通用执行单元,也叫 CUDA core,对应逻辑结构的一个 thread,它负责执行具体的加减乘除、加载存储等运算指令。

Tensor Core:专用执行单元,专用于矩阵乘法

整个 SM 内部的协作过程:Instruction Cache 提供已经缓存的指令,减少从全局内存加载指令的延迟,确保 Warp Scheduler能够快速获取所需指令。Warp Scheduler 从 Instruction Cache 中提取指令,根据当前的warp状态和依赖关系,选择合适的warp进行执行。Dispatch Unit 从Warp Scheduler那里获取指令,并将其分发到 执行单元。Instruction Cache的作用是让Dispatch Unit能够快速获得已经缓存的指令,减少指令调度延迟。

5.GPU 内存层次及每种存储的特点

首先有个计算机基本概念,距离执行计算单元越近的内存越宝贵,读写速度越快,尺寸也越小。

  • 全局内存(Global Memory):GPU 上最大的内存,即常说的 HBM 内存,可以被所有线程访问,从 host 端搬运的数据存放在此处。
  • L2 Cache:SM 访问交互的缓存存放区域,即会放 global/local memory 的数据,当 thread 需要数据时,会先访问 L2 Cache,若没有才访问 global memory, 这时候的读取不会只读单独的数据,而是把周围一部分数据一起带到 L2 Cache上,方便之后或者临近的 thread 能够在L2 Cache上找到它(局部性原理)。
  • 共享内存(Shared Memory)/L1 Cache:位于 sm 内部,同一个 block 的线程可访问同一块共享内存,在物理结构上和L1 cache 位于同一块物理内存,但 shmem 可由程序员手动管理, L1 cache 由硬件自动管理。
  • 本地内存(Local Memory):每个线程都可以使用自己的本地内存,可以在其中存储临时变量。它具有最小的范围并且专用于每个单独的线程。
  • Register Cache:寄存器是每个线程私有的,访问速度最快,每个SM中,寄存器缓存(Register Cache)用于提供线程的寄存器访问。每个SM包含一定数量的寄存器(physical registers),所有线程都可以从这些寄存器中分配资源。

6.CUDA stream

stream 相等于是 GPU 上的任务队列,即命令流水线,它允许多个 CUDA 操作在不同的 stream 中并行执行,从而提高 GPU 的利用率和性能。一个 stream 上的任务是严格顺序执行的,上一个命令执行结束才会执行下一个命令,不同 stream 的命令之间无执行顺序,即可以并发执行不同命令。

7.CPU、GPU、NPU的区别

CPU 有更多的控制单元,更大的缓存,更少的 ALU(算术逻辑单元),延迟较低;相较之下,GPU 使用大量的 ALU 来最大化计算能力和吞吐量,少量的缓存和控制单元,延迟较高。

所以对于他们的设计目标来说,CPU 被设计为顺序执行指令,为了提高顺序执行性能,多年来 CPU 设计中引入了许多功能,包括指令流水线、乱序执行、推测执行和多级缓存等。CPU 的重点是减少指令执行延迟,以便 CPU 能够尽快执行指令序列。

GPU 专为大规模并行性和高吞吐量而设计,但代价是较高的指令延迟。

NPU,即神经网络加速器,在电路层模拟神经元,相比于 CPU 和 GPU,NPU 通过突出权重实现存储和计算一体化,从而提高运行效率,一条指令完成一组神经元的处理,提高运行效率

8.DMA RDMA

DMA(直接内存访问):传统的内存访问需要通过 CPU 进行数据 copy 来移动数据,通过 CPU 将内存中的 buffer1 移动到 buffer2,而 DMA 模式在 DMA Engine 之间通过硬件直接将 buffer1 移动到 buffer2,不需要 cpu 的参与,减少了 copy 开销。

RDMA:即远程直接内存访问,指多个计算机之间访问使用 DMA,一个主机可以直接访问另一个主机的内存。

使用 RDMA 的优势如下:

  • 零拷贝(Zero-copy) - 应用程序能够直接执行数据传输,在不涉及到网络软件栈的情况下。数据能够被直接发送到缓冲区或者能够直接从缓冲区里接收,而不需要被复制到网络层。
  • 内核旁路(Kernel bypass) - 应用程序可以直接在用户态执行数据传输,不需要在内核态与用户态之间做上下文切换。
  • 不需要 CPU 干预(No CPU involvement) - 应用程序可以访问远程主机内存而不消耗远程主机中的任何 CPU。远程主机内存能够被读取而不需要远程主机上的进程(或 CPU)参与。远程主机的 CPU 的缓存(cache)不会被访问的内存内容所填充。
  • 消息基于事务(Message based transactions) - 数据被处理为离散消息而不是流,消除了应用程序将流切割为不同消息/事务的需求。
  • 支持分散/聚合条目(Scatter/gather entries support) - RDMA 原生态支持分散/聚合。也就是说,读取多个内存缓冲区然后作为一个流发出去或者接收一个流然后写入到多个内存缓冲区里去。

waiting…