I have written a simple CUDA program to perform array reduction using thread block clusters and distributed shared memory. I am compiling it with CUDA 12.0 and running on a hopper GPU. Below is the code I use:
#include <stdio.h>
#include <cooperative_groups.h>
#define CLUSTER_SIZE 4
#define BLOCK_SIZE 32
namespace cg = cooperative_groups;
__global__ void __cluster_dims__(CLUSTER_SIZE, 1, 1)
cluster_reduce_sum(int n, float *arr, float *sum)
{
__shared__ float shared_mem[BLOCK_SIZE];
__shared__ float cluster_sum;
cluster_sum = 0.0f;
cg::cluster_group cluster = cg::this_cluster();
unsigned int cluster_block_rank = cluster.block_rank();
unsigned int cluster_size = cluster.dim_blocks().x;
const int idx = blockIdx.x * blockDim.x + threadIdx.x;
shared_mem[threadIdx.x] = 0.0f;
if (idx < n) {
shared_mem[threadIdx.x] = arr[idx];
}
__syncthreads();
for (int offset = BLOCK_SIZE / 2; offset; offset /= 2) {
if (threadIdx.x < offset) {
shared_mem[threadIdx.x] += shared_mem[threadIdx.x + offset];
}
__syncthreads();
}
cluster.sync();
if (threadIdx.x == 0) {
atomicAdd(cluster.map_shared_rank(&cluster_sum, 0), shared_mem[0]);
printf("blockIdx: %d cluster_sum: %f\n", blockIdx.x, (float)*cluster.map_shared_rank(&cluster_sum, 0));
}
cluster.sync();
if (threadIdx.x == 0 && cluster_block_rank == 0) {
atomicAdd(sum, cluster_sum);
}
cluster.sync();
}
int main(int argc, char* argv[]) {
int n = 128;
if (argc > 1) {
n = atoi(argv[1]);
}
float *h_arr, *h_sum, sum;
h_arr = (float*) malloc(n * sizeof(float));
h_sum = (float*) malloc(sizeof(float));
int upper = 1024, lower = -1024;
sum = 0.0f;
for(int i = 0; i < n; i++)
{
h_arr[i] = 1;
sum += h_arr[i];
}
float *d_arr, *d_sum;
cudaMalloc(&d_arr, n * sizeof(float));
cudaMalloc(&d_sum, sizeof(float));
cudaMemcpy(d_arr, h_arr, n * sizeof(float), cudaMemcpyHostToDevice);
cudaMemset(d_sum, 0, sizeof(float));
int num_clusters = (n-1) / (CLUSTER_SIZE * BLOCK_SIZE) + 1;
cluster_reduce_sum <<< CLUSTER_SIZE * num_clusters, BLOCK_SIZE >>> (n, d_arr, d_sum);
cudaError_t error = cudaGetLastError();
if(error != cudaSuccess)
{
printf("CUDA error: %s\n", cudaGetErrorString(error));
return -1;
}
cudaMemcpy(h_sum, d_sum, sizeof(float), cudaMemcpyDeviceToHost);
if (*h_sum != sum) {
printf("Kernel incorrect: %f vs %f\n", sum, *h_sum);
}
return 0;
}
The outputs of the kernel and CPU do not match and as far as I know, there doesn’t seem to be any bugs in my code. Please help me find the issue with this code. Thanks!
Edit: including some debug info to the question. I added a print statement just before adding the per block sum (shared_mem[0]
)to the cluster_sum
. The per block sum seem to be all correct (32). However when I print the final cluster just before adding it to the total sum, it is not correct. Expected value of cluster_sum
is 128, but I see 64. Seems the block sum is not properly accumulated to the cluster sum.
Another interesting point I observed is that if I just change the input datatype to int
from float
, the code works completely fine and passes the output check.
Edit: I run the code. It looks that atomicAdd does not work correctly on float. I add several printf to the original codes.
atomicAdd
. Probably unrelated, but what is the purpose of the lastcluster.sync()
?shared_mem[0]
just from the first thread of the first block of the cluster (i.e. gather instead of scatter the block results). That way you don't need distributed shared memory atomics (Are they supported? No mention of them in the atomics chapter of the programming guide, I think).