LeetGPU | Reduction 为例,这是一个典型的 Memory-Bound 的操作。

注:本文的性能测试数据为 NVIDIA GeForce RTX 4060 Laptop 运行结果,仅供参考。

单线程串行求和

只有一个 GPU 线程工作

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#include <cuda_runtime.h>

__global__ void reduceKernel(const float *input, float *output, int N) {
double sum = 0.0;
for (int i = 0; i < N; ++i) {
sum += input[i];
}

output[0] = static_cast<float>(sum);
}

extern "C" void solve(const float *input, float *output, int N) {
cudaMemset(output, 0, sizeof(float));
if (N <= 0) {
return;
}

reduceKernel<<<1, 1>>>(input, output, N);
}

CPU 将 reduceKernel 交给 GPU,然后立即返回执行后面的代码;GPU 收到任务后慢慢执行,利用 nsys 工具观察 CUDA API 的时间线,会发现 CPU 侧调用 CUDA Runtime/Driver,将 kernel launch 命令提交出去的时间仅为 60.145us,接下来 cudaDeviceSynchronize() 约为 105.407ms,相当于 CPU 在等待 GPU 完成任务。

多线程并行读取

最直接的并行化方式是每个线程处理若干个元素,每个线程计算一个局部和,然后使用 atomicAdd 加到最终结果。

相比单线程版本,它能够并行读取输入,但是每个线程都会进行一次原子操作,竞争严重。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#include <cuda_runtime.h>

__global__ void reduceKernel(const float *input, float *output, int N) {
int idx = blockDim.x * blockIdx.x + threadIdx.x;
float sum = 0.0f;
sum += input[idx];
atomicAdd(output, sum);
}

extern "C" void solve(const float *input, float *output, int N) {
cudaMemset(output, 0, sizeof(float));
if (N <= 0) {
return;
}

int ThreadsPerBlock = 256;
int BlocksPerGrid = (N + ThreadsPerBlock - 1) / ThreadsPerBlock;

reduceKernel<<<BlocksPerGrid, ThreadsPerBlock>>>(input, output, N);
}

注:精度错误。

块内规约

朴素树形

最直接的块内树形规约是从步长 1 开始,每一轮让满足条件的线程把右侧元素累加到当前位置。每个线程将局部计算的结果写入共享内存,每个 Block 内先完成树形规约,这样就可以将原子加操作数量从线程数下降到块数。不过规约的每一层都需要执行 __syncthreads()

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
#include <cuda_runtime.h>

#define BLOCK_SIZE 256

__global__ void reduceKernel(const float *input, float *output, int N) {
int idx = blockDim.x * blockIdx.x + threadIdx.x;
int tid = threadIdx.x;
__shared__ float shared[BLOCK_SIZE];

shared[tid] = idx < N ? input[idx] : 0.0f;
__syncthreads();

for (int stride = 1; stride < BLOCK_SIZE; stride <<= 1) {
if (tid % (2 * stride) == 0) {
shared[tid] += shared[tid + stride];
}
__syncthreads();
}

if (tid == 0) {
atomicAdd(output, shared[0]);
}
}

extern "C" void solve(const float *input, float *output, int N) {
cudaMemset(output, 0, sizeof(float));
if (N <= 0) {
return;
}

int ThreadsPerBlock = BLOCK_SIZE;
int BlocksPerGrid = (N + ThreadsPerBlock - 1) / ThreadsPerBlock;

reduceKernel<<<BlocksPerGrid, ThreadsPerBlock>>>(input, output, N);
}

随着 stride 增大,同一个 Warp 内满足 threadIdx.x % (2 * stride) == 0 的线程越来越稀疏,会产生明显的 Warp Divergence;同时取模运算和交错的 Shared Memory 访问也会带来额外开销。

例如

  • stride=1 时:每个 Warp 内只有偶数线程工作 → 50% 利用率
  • stride=2 时:每 4 个线程只有 1 个工作 → 25% 利用率
  • ……

在 ncu 的 Memory Chart 中

1
2
Grid Size  = 16,384 blocks
Block Size = 256 threads
  • Global

input[idx] 是 global read,每个 warp 执行一次 global load instruction,因此总数为 atomicAdd(output, shared[0]); 每个 block 执行一次,共 ,总数即为 .

  • Shared Store

开始写入 shared,每个 block 8 个 warp

后续规约的循环中,随着 stride 递增,活跃的 warp 总数为 ,即每个 block 的 reduction 需要由 47 个 warp-instance 执行。整个 grid 总共 .

  • Shared Load

shared[tid] += shared[tid + stride]; 每个活跃的 warp 两次 shared load,最后原子加还需要一次,所以一个 block 共 次 load,整个 grid 共

  • Device Memory

总线程数为 ,每个线程读取 float,总数即为 .

交错寻址

可以把线程是否参与规约的判断改成连续线程计算目标下标。前半部分线程通过 index = 2 * stride * tid 找到本轮需要合并的位置,参与计算的线程集中在 Warp 前部,减少了朴素树形中由取模条件造成的 Warp Divergence。

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
#include <cuda_runtime.h>

#define BLOCK_SIZE 256

__global__ void reduceKernel(const float *input, float *output, int N) {
int tid = threadIdx.x;
int idx = blockDim.x * blockIdx.x + tid;
__shared__ float shared[BLOCK_SIZE];

shared[tid] = idx < N ? input[idx] : 0.0f;
__syncthreads();

for (int stride = 1; stride < BLOCK_SIZE; stride <<= 1) {
int index = 2 * stride * tid;
if (index < BLOCK_SIZE) {
shared[index] += shared[index + stride];
}
__syncthreads();
}

if (tid == 0) {
atomicAdd(output, shared[0]);
}
}

extern "C" void solve(const float *input, float *output, int N) {
cudaMemset(output, 0, sizeof(float));
if (N <= 0) {
return;
}

int BlocksPerGrid = (N + BLOCK_SIZE - 1) / BLOCK_SIZE;
reduceKernel<<<BlocksPerGrid, BLOCK_SIZE>>>(input, output, N);
}

不过随着 stride 增大,Shared Memory 的访问间隔也同步增大,多个线程可能映射到相同的 Memory Bank,因此仍然可能出现 Bank Conflict。

例如

  • 第 1 轮循环(stride = 1):tid 0 访问 shared [0]、shared [1];tid 16 访问 shared [32]、shared [33]。shared [0] 和 shared [32] 都落在 Bank 0 → 2 路 Bank Conflict
  • 第 2 轮循环(stride = 2):tid 0 访问 shared [0,2];tid 8 访问 shared [32,34];tid 16 访问 shared [64,66];tid 24 访问 shared [96,98]。shared [0]、shared [32]、shared [64]、shared [96] 都在 Bank 0 → 4 路 Bank Conflict
  • 第 3 轮循环(stride = 4):8 路 Bank Conflict
  • ……

步长反转

把规约方向反过来,从 BLOCK_SIZE / 2 开始不断将步长减半。每一轮都由连续的前 offset 个线程访问 shared[tid]shared[tid + offset],避免了交错寻址中的跨步 Shared Memory 访问,也让有效线程连续分布,从而消除典型的 Bank Conflict 和主要的 Warp Divergence。最后一个 Warp 中活跃线程仍会逐步减少,后续再通过展开最后 Warp 继续优化。

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
#include <cuda_runtime.h>

#define BLOCK_SIZE 256

__global__ void reduceKernel(const float *input, float *output, int N) {
int idx = blockDim.x * blockIdx.x + threadIdx.x;
int tid = threadIdx.x;
__shared__ float shared[BLOCK_SIZE];

shared[tid] = idx < N ? input[idx] : 0.0f;
__syncthreads();

for (int stride = BLOCK_SIZE / 2; stride > 0; stride >>= 1) {
if (tid < stride) {
shared[tid] += shared[tid + stride];
}
__syncthreads();
}

if (tid == 0) {
atomicAdd(output, shared[0]);
}
}

extern "C" void solve(const float *input, float *output, int N) {
cudaMemset(output, 0, sizeof(float));
if (N <= 0) {
return;
}

int ThreadsPerBlock = BLOCK_SIZE;
int BlocksPerGrid = (N + ThreadsPerBlock - 1) / ThreadsPerBlock;

reduceKernel<<<BlocksPerGrid, ThreadsPerBlock>>>(input, output, N);
}
  • 活跃线程 → Warp Divergence
stride 朴素树形 交错寻址 stride 步长反转
1 0, 2, 4, ..., 254 0, 1, 2, ..., 127 128 0, 1, 2, ..., 127
2 0, 4, 8, ..., 252 0, 1, 2, ..., 63 64 0, 1, 2, ..., 63
3 0, 8, 16, ..., 248 0, 1, 2, ..., 31 32 0, 1, 2, ..., 31
4 0, 16, 32, ..., 240 0, 1, 2, ..., 15 16 0, 1, 2, ..., 15
  • 同一个 warp 访问共享内存地址与 Bank → Bank Conflict
规约轮次 朴素树形 交错寻址 步长反转 Bank Conflict
1 (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), (30,31) (0,1), (2,3),..., (30,31) // (32,33), (34,35),..., (62,63) (0,128), (1,129), (2,130),..., (31,159) 交错寻址 2-way;朴素树形、步长反转无冲突
2 (0,2), (4,6), (8,10), (12,14), (16,18), (20,22), (24,26), (28,30) (0,2), (4,6),..., (28,30) // (32,34), (36,38),..., (60,62) // (64,66), (68,70),..., (92,94) // (96,98), (100,102),..., (124,126) (0,64), (1,65), (2,66),..., (31,95) 交错寻址 4-way;朴素树形、步长反转无冲突
3 (0,4), (8,12), (16,20), (24, 28) (0,4), (8,12), (16,20), (24,28) // (32,36), (40,44), (48,52), (56,60) // (64,68), (72,76), (80,84), (88,92) // (96,100), (104,108), (112,116), (120,124) // (128,132), (136,140), (144,148), (152,156) // (160,164), (168,172), (176,180), (184,188) // (192,196), (200,204), (208,212), (216,220) // (224,228), (232,236), (240,244), (248,252) (0,32), (1,33), (2,34),..., (31,63) 交错寻址 8-way;朴素树形、步长反转无冲突
4 (0,8), (16,24) (0,8), (16,24) // (32,40), (48,56) // (64,72), (80,88) // (96,104), (112,120) // (128,136), (144,152) // (160,168), (176,184) // (192,200), (208,216) // (224,232), (240,248) (0,16), (1,17), (2,18),..., (15,31) 交错寻址 16-way;朴素树形、步长反转无冲突
5 (0,16) (0,16) // (32,48) // (64,80) // (96,112) // (128,144) // (160,176) // (192,208) // (224,240) (0,8), (1,9), (2,10), (3,11), (4,12), (5,13), (6,14), (7,15) 交错寻址最高可形成 8-way 有效冲突;朴素树形、步长反转无冲突

双元素处理

前面的版本每个线程只在进入 Shared Memory 规约前读取一个元素。可以让每个线程先在寄存器中累加两个输入元素,再把局部和写入 Shared Memory。这样同样数量的线程一次可以覆盖两倍的数据,Block 数量近似减半,同时前几轮本来会逐步退出的线程被提前用于 Global Memory 读取和加法。

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
#include <cuda_runtime.h>

#define BLOCK_SIZE 256

__global__ void reduceKernel(const float *input, float *output, int N) {
int idx = 2 * blockDim.x * blockIdx.x + threadIdx.x;
int tid = threadIdx.x;
__shared__ float shared[BLOCK_SIZE];

double sum = 0.0f;
if (idx < N) {
sum += input[idx];
}
if (idx + BLOCK_SIZE < N) {
sum += input[idx + BLOCK_SIZE];
}

shared[tid] = sum;
__syncthreads();

for (int stride = BLOCK_SIZE / 2; stride > 0; stride >>= 1) {
if (tid < stride) {
shared[tid] += shared[tid + stride];
}
__syncthreads();
}

if (tid == 0) {
atomicAdd(output, shared[0]);
}
}

extern "C" void solve(const float *input, float *output, int N) {
cudaMemset(output, 0, sizeof(float));
if (N <= 0) {
return;
}

int ThreadsPerBlock = BLOCK_SIZE;
int elementsPerBlock = BLOCK_SIZE * 2;
int BlocksPerGrid = (N + elementsPerBlock - 1) / elementsPerBlock;

reduceKernel<<<BlocksPerGrid, ThreadsPerBlock>>>(input, output, N);
}

展开最后 Warp

当规约只剩最后 32 个线程时,后续计算全部发生在同一个 Warp 内。此时不再需要每一层都执行 Block 级别的 __syncthreads(),可以显式展开最后一个 Warp 的加法,只保留 Warp 级同步。

在 Volta 及之后的架构中存在 Independent Thread Scheduling,因此下面使用 __syncwarp() 保证同一个 Warp 对 Shared Memory 的读写顺序,而不是依赖早期架构隐式的 Warp Lockstep 行为。

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
#include <cuda_runtime.h>

#define BLOCK_SIZE 256

__global__ void reduceKernel(const float *input, float *output, int N) {
int idx = 2 * blockDim.x * blockIdx.x + threadIdx.x;
int tid = threadIdx.x;
__shared__ float shared[BLOCK_SIZE];

double sum = 0.0f;
if (idx < N) {
sum += input[idx];
}
if (idx + BLOCK_SIZE < N) {
sum += input[idx + BLOCK_SIZE];
}

shared[tid] = sum;
__syncthreads();

for (int stride = BLOCK_SIZE / 2; stride > 32; stride >>= 1) {
if (tid < stride) {
shared[tid] += shared[tid + stride];
}
__syncthreads();
}

if (tid < 32) {
volatile float *v = shared;
v[tid] += v[tid + 32];
__syncwarp();
v[tid] += v[tid + 16];
__syncwarp();
v[tid] += v[tid + 8];
__syncwarp();
v[tid] += v[tid + 4];
__syncwarp();
v[tid] += v[tid + 2];
__syncwarp();
v[tid] += v[tid + 1];
}

if (tid == 0) {
atomicAdd(output, shared[0]);
}
}

extern "C" void solve(const float *input, float *output, int N) {
cudaMemset(output, 0, sizeof(float));
if (N <= 0) {
return;
}

int ThreadsPerBlock = BLOCK_SIZE;
int elementsPerBlock = BLOCK_SIZE * 2;
int BlocksPerGrid = (N + elementsPerBlock - 1) / elementsPerBlock;

reduceKernel<<<BlocksPerGrid, ThreadsPerBlock>>>(input, output, N);
}

完全循环展开

BLOCK_SIZE 在 Kernel 启动时通常是固定值,因此可以把它变成模板参数。编译器可以在编译期确定每一层规约是否存在,并把 Shared Memory 规约和最后一个 Warp 的规约全部展开,消除循环变量、循环条件和分支控制带来的指令开销。

下面每个 if 的条件都只依赖模板参数,是编译期常量。对于 reduceKernel<256>,编译器会进行常量折叠,只保留 256 线程版本实际需要的规约步骤,因此源码中不再存在运行时规约循环。

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
#include <cuda_runtime.h>

#define BLOCK_SIZE 256

template <unsigned int BLOCK_SIZE_>
__global__ void reduceKernel(const float *input, float *output, int N) {
int idx = 2 * BLOCK_SIZE_ * blockIdx.x + threadIdx.x;
int tid = threadIdx.x;
__shared__ float shared[BLOCK_SIZE_];

double sum = 0.0f;
if (idx < N) {
sum += input[idx];
}
if (idx + BLOCK_SIZE_ < N) {
sum += input[idx + BLOCK_SIZE_];
}

shared[tid] = sum;
__syncthreads();

if (BLOCK_SIZE_ >= 1024) {
if (tid < 512)
shared[tid] += shared[tid + 512];
__syncthreads();
}
if (BLOCK_SIZE_ >= 512) {
if (tid < 256)
shared[tid] += shared[tid + 256];
__syncthreads();
}
if (BLOCK_SIZE_ >= 256) {
if (tid < 128)
shared[tid] += shared[tid + 128];
__syncthreads();
}
if (BLOCK_SIZE_ >= 128) {
if (tid < 64)
shared[tid] += shared[tid + 64];
__syncthreads();
}

if (tid < 32) {
volatile float *v = shared;
v[tid] += v[tid + 32];
__syncwarp();
v[tid] += v[tid + 16];
__syncwarp();
v[tid] += v[tid + 8];
__syncwarp();
v[tid] += v[tid + 4];
__syncwarp();
v[tid] += v[tid + 2];
__syncwarp();
v[tid] += v[tid + 1];
}

if (tid == 0) {
atomicAdd(output, shared[0]);
}
}

extern "C" void solve(const float *input, float *output, int N) {
cudaMemset(output, 0, sizeof(float));
if (N <= 0) {
return;
}

int ThreadsPerBlock = BLOCK_SIZE;
int elementsPerBlock = BLOCK_SIZE * 2;
int BlocksPerGrid = (N + elementsPerBlock - 1) / elementsPerBlock;

reduceKernel<BLOCK_SIZE>
<<<BlocksPerGrid, ThreadsPerBlock>>>(input, output, N);
}

Warp Shuffle

Shared Memory 主要用于不同线程之间交换中间结果,而同一个 Warp 内可以直接通过 Shuffle 指令交换寄存器值。使用 __shfl_down_sync 后,每个 Warp 的规约都在寄存器中完成,只需要让每个 Warp 的 Lane 0 把局部和写入一小段 Shared Memory,再由第一个 Warp 完成 Block 级规约。

这样 Shared Memory 从 BLOCK_SIZE 个元素缩小到 BLOCK_SIZE / 32 个元素,并且整个 Block 只需要一次 __syncthreads()

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
#include <cuda_runtime.h>

#define BLOCK_SIZE 256

__device__ __forceinline__ float warpReduceSum(float value) {
value += __shfl_down_sync(0xffffffff, value, 16);
value += __shfl_down_sync(0xffffffff, value, 8);
value += __shfl_down_sync(0xffffffff, value, 4);
value += __shfl_down_sync(0xffffffff, value, 2);
value += __shfl_down_sync(0xffffffff, value, 1);
return value;
}

__global__ void reduceKernel(const float *input, float *output, int N) {
int idx = blockDim.x * blockIdx.x * 2 + threadIdx.x;
int tid = threadIdx.x;
__shared__ float shared[BLOCK_SIZE / 32];

double sum = 0.0f;
if (idx < N) {
sum += input[idx];
}
if (idx + BLOCK_SIZE < N) {
sum += input[idx + BLOCK_SIZE];
}

sum = warpReduceSum(sum);

int lane_id = tid % 32;
int warp_id = tid >> 5;

if (lane_id == 0) {
shared[warp_id] = sum;
}
__syncthreads();

if (warp_id == 0) {
double warp_sum = lane_id < BLOCK_SIZE / 32 ? shared[lane_id] : 0.0f;
double block_sum = warpReduceSum(warp_sum);
if (lane_id == 0) {
atomicAdd(output, block_sum);
}
}
}

extern "C" void solve(const float *input, float *output, int N) {
cudaMemset(output, 0, sizeof(float));
if (N <= 0) {
return;
}

int ThreadsPerBlock = BLOCK_SIZE;
int elementsPerBlock = ThreadsPerBlock * 2;
int BlocksPerGrid = (N + elementsPerBlock - 1) / elementsPerBlock;

reduceKernel<<<BlocksPerGrid, ThreadsPerBlock>>>(input, output, N);
}

网格步长循环

每个线程不再只读取一两个元素,而是按整个网格的跨度连续处理多个元素。这样可以使用固定数量的块处理任意长度的输入,每个线程循环处理多个元素。

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
#include <cuda_runtime.h>

#define BLOCK_SIZE 256

__device__ __forceinline__ float warpReduceSum(float value) {
value += __shfl_down_sync(0xffffffff, value, 16);
value += __shfl_down_sync(0xffffffff, value, 8);
value += __shfl_down_sync(0xffffffff, value, 4);
value += __shfl_down_sync(0xffffffff, value, 2);
value += __shfl_down_sync(0xffffffff, value, 1);
return value;
}

__global__ void reduceKernel(const float *input, float *output, int N) {
int idx = blockDim.x * blockIdx.x * 2 + threadIdx.x;
int stride = gridDim.x * blockDim.x * 2;
int tid = threadIdx.x;
__shared__ float shared[BLOCK_SIZE / 32];

double sum = 0.0f;
for (int i = idx; i < N; i += stride) {
sum += input[i];
if (i + BLOCK_SIZE < N) {
sum += input[i + BLOCK_SIZE];
}
}
sum = warpReduceSum(sum);

int lane_id = tid % 32;
int warp_id = tid >> 5;

if (lane_id == 0) {
shared[warp_id] = sum;
}
__syncthreads();

if (warp_id == 0) {
double warp_sum = lane_id < BLOCK_SIZE / 32 ? shared[lane_id] : 0.0f;
double block_sum = warpReduceSum(warp_sum);
if (lane_id == 0) {
atomicAdd(output, block_sum);
}
}
}

extern "C" void solve(const float *input, float *output, int N) {
cudaMemset(output, 0, sizeof(float));
if (N <= 0) {
return;
}

int ThreadsPerBlock = BLOCK_SIZE;
int elementsPerBlock = ThreadsPerBlock * 2;
int BlocksPerGrid = (N + elementsPerBlock - 1) / elementsPerBlock;
BlocksPerGrid = BlocksPerGrid > 1024 ? 1024 : BlocksPerGrid;

reduceKernel<<<BlocksPerGrid, ThreadsPerBlock>>>(input, output, N);
}

向量化 + Stride Loop

Grid-Stride Loop 解决的是数据规模大于当前 Grid 覆盖范围的问题,而向量化读取可以减少 Global Memory Load 指令数量并提高单线程每次访存的有效载荷。将连续的 4 个 float 视为一个 float4 后,每个线程一次读取 16 Byte,并继续按照整个 Grid 的跨度循环读取后续向量。

这种方式既保持了连续线程之间的合并访存,也允许固定数量的 Block 反复处理任意长度的输入。向量部分结束后再单独处理不足 4 个元素的尾部,因此不会遗漏数据。这里默认 input 具有至少 16 Byte 对齐,cudaMalloc 返回的设备地址满足这一要求。

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
#include <cuda_runtime.h>

#define BLOCK_SIZE 256

__device__ __forceinline__ float warpReduceSum(float value) {
value += __shfl_down_sync(0xffffffff, value, 16);
value += __shfl_down_sync(0xffffffff, value, 8);
value += __shfl_down_sync(0xffffffff, value, 4);
value += __shfl_down_sync(0xffffffff, value, 2);
value += __shfl_down_sync(0xffffffff, value, 1);
return value;
}

__global__ void reduceKernel(const float *input, float *output, int N) {
int idx = blockDim.x * blockIdx.x + threadIdx.x;
int stride = gridDim.x * blockDim.x;
int tid = threadIdx.x;
__shared__ float shared[BLOCK_SIZE / 32];

double sum = 0.0f;
const float4 *input4 = reinterpret_cast<const float4 *>(input);
for (int i = idx; i < N / 4; i += stride) {
float4 value = input4[i];
sum += value.x + value.y + value.z + value.w;
}

int tail = N / 4 * 4;
for (int i = idx + tail; i < N; i += stride) {
sum += input[i];
}

sum = warpReduceSum(sum);

int lane_id = tid % 32;
int warp_id = tid >> 5;

if (lane_id == 0) {
shared[warp_id] = sum;
}
__syncthreads();

if (warp_id == 0) {
double warp_sum = lane_id < BLOCK_SIZE / 32 ? shared[lane_id] : 0.0f;
double block_sum = warpReduceSum(warp_sum);
if (lane_id == 0) {
atomicAdd(output, block_sum);
}
}
}

extern "C" void solve(const float *input, float *output, int N) {
cudaMemset(output, 0, sizeof(float));
if (N <= 0) {
return;
}

int ThreadsPerBlock = BLOCK_SIZE;
int elementsPerBlock = ThreadsPerBlock;
int BlocksPerGrid = (N + elementsPerBlock - 1) / elementsPerBlock;
BlocksPerGrid = BlocksPerGrid > 1024 ? 1024 : BlocksPerGrid;

reduceKernel<<<BlocksPerGrid, ThreadsPerBlock>>>(input, output, N);
}