Illegal memory access in reduce kernel when built with CUDA Toolkit 13.1 #3109

Activity

1. zcbenz

   added

   bug

   cuda

   on Feb 7, 2026
2. zcbenz commented on Feb 10, 2026

   

   zcbenz

   on Feb 10, 2026

   Last edited by zcbenz

   MemberAuthor

   More actions

   Independent test below. Note that the CUDA arch must be set to 120a.

   Details

   // bug.cu
   // Build: nvcc bug.cu -o bug -std=c++20 "--generate-code=arch=compute_120a,code=[compute_120a,sm_120a]"
   // Run: ./bug
   
   #include <cuda_runtime.h>
   #include <iostream>
   #include <cstdlib>
   #include <cuda/std/array>
   
   #include <cooperative_groups.h>
   #include <cooperative_groups/reduce.h>
   #include <cub/block/block_load.cuh>
   #include <cub/cub.cuh>
   
   namespace {
   
   namespace cg = cooperative_groups;
   
   #define WARP_SIZE 32
   
   struct ColReduceArgs {
     size_t reduction_size;
     int64_t reduction_stride;
   
     cuda::std::array<int32_t, 10> reduce_shape;
     cuda::std::array<int64_t, 10> reduce_strides;
     int reduce_ndim;
   
     size_t non_col_reductions;
   };
   
   template <typename IdxT = int64_t>
   inline __host__ __device__ IdxT
   elem_to_loc(IdxT elem, const int* shape, const int64_t* strides, int ndim) {
     IdxT loc = 0;
     for (int i = ndim - 1; i >= 0 && elem > 0; --i) {
       loc += (elem % shape[i]) * IdxT(strides[i]);
       elem /= shape[i];
     }
     return loc;
   }
   
   template <int DIM, bool General = true, typename OffsetT = size_t>
   struct LoopedElemToLoc {
     int dim;
     LoopedElemToLoc<DIM - 1, General, OffsetT> inner_looper;
     OffsetT offset{0};
     int index{0};
   
     __device__ LoopedElemToLoc(int dim) : dim(dim), inner_looper(dim - 1) {}
   
     __device__ void next(const int* shape, const int64_t* strides) {
       if (dim == 0) {
         return;
       }
       index++;
       offset += OffsetT(strides[dim - 1]);
       if (index >= shape[dim - 1]) {
         index = 0;
         inner_looper.next(shape, strides);
         offset = inner_looper.offset;
       }
     }
   
     __device__ void next(int n, const int* shape, const int64_t* strides) {
       if (dim == 0) {
         return;
       }
       index += n;
       offset += n * OffsetT(strides[dim - 1]);
   
       if (index >= shape[dim - 1]) {
         int extra = index - shape[dim - 1];
         if (extra >= shape[dim - 1]) {
           inner_looper.next(1 + extra / shape[dim - 1], shape, strides);
           extra = extra % shape[dim - 1];
         } else {
           inner_looper.next(shape, strides);
         }
         index = 0;
         offset = inner_looper.offset;
         if (extra > 0) {
           next(extra, shape, strides);
         }
       }
     }
   
     __device__ OffsetT location() {
       return offset;
     }
   };
   
   template <typename OffsetT>
   struct LoopedElemToLoc<1, true, OffsetT> {
     int dim;
     OffsetT offset{0};
     int index{0};
   
     __device__ LoopedElemToLoc(int dim) : dim(dim) {}
   
     __device__ void next(const int* shape, const int64_t* strides) {
       index++;
       if (dim > 1) {
         offset = elem_to_loc<OffsetT>(index, shape, strides, dim);
       } else {
         offset += OffsetT(strides[0]);
       }
     }
   
     __device__ void next(int n, const int* shape, const int64_t* strides) {
       index += n;
       if (dim > 1) {
         offset = elem_to_loc<OffsetT>(index, shape, strides, dim);
       } else {
         offset = index * OffsetT(strides[0]);
       }
     }
   
     __device__ OffsetT location() {
       return offset;
     }
   };
   
   template <typename OffsetT>
   struct LoopedElemToLoc<1, false, OffsetT> {
     OffsetT offset{0};
   
     __device__ LoopedElemToLoc(int) {}
   
     __device__ void next(const int*, const int64_t* strides) {
       offset += OffsetT(strides[0]);
     }
   
     __device__ void next(int n, const int*, const int64_t* strides) {
       offset += n * OffsetT(strides[0]);
     }
   
     __device__ OffsetT location() {
       return offset;
     }
   };
   
   void checkCuda(cudaError_t err, const char* where) {
     if (err != cudaSuccess) {
       std::cerr << "CUDA error at " << where << ": " << cudaGetErrorString(err) << "\n";
       std::exit(1);
     }
   }
   
   }  // namespace
   
   template <
       typename T,
       int NDIM,
       int BM,
       int BN,
       int N_READS = 4>
   __global__ void col_reduce_looped(
       T* in,
       const __grid_constant__ ColReduceArgs args) {
     auto grid = cg::this_grid();
     auto block = cg::this_thread_block();
     auto warp = cg::tiled_partition<WARP_SIZE>(block);
   
     constexpr int threads_per_row = BN / N_READS;
   
     size_t tile_idx = grid.block_rank();
     size_t tile_x = tile_idx % ((args.reduction_stride + BN - 1) / BN);
   
     short thread_x = block.thread_rank() % threads_per_row;
     short thread_y = block.thread_rank() / threads_per_row;
   
     T totals[N_READS];
     for (int i = 0; i < N_READS; i++) {
       totals[i] = 0;
     }
   
     size_t total = args.non_col_reductions * args.reduction_size;
     LoopedElemToLoc<NDIM, (NDIM > 2)> loop(args.reduce_ndim);
     loop.next(thread_y, args.reduce_shape.data(), args.reduce_strides.data());
     for (size_t r = thread_y; r < total; r += BM) {
       T vals[N_READS];
       cub::LoadDirectBlocked(
           thread_x,
           in + loop.location(),
           vals,
           args.reduction_stride - tile_x * BN,
           0);
       for (int i = 0; i < N_READS; i++) {
         totals[i] = vals[i];
       }
       loop.next(BM, args.reduce_shape.data(), args.reduce_strides.data());
     }
   
     constexpr int n_outputs = BN / threads_per_row;
     __shared__ T shared_vals[BM * BN];
     short s_idx = thread_y * BN + thread_x * N_READS;
     for (int i = 0; i < N_READS; i++) {
       shared_vals[s_idx + i] = totals[i];
     }
     block.sync();
     for (int i = 0; i < n_outputs; i++) {
       totals[i] = shared_vals[s_idx + i];
     }
   }
   
   int main() {
     float* a;
     checkCuda(cudaMalloc(&a, 16 * sizeof(float)), "cudaMallocManaged");
   
     ColReduceArgs args{
         2, 2,
         {2, 2},
         {8, 2},
         2,
         2};
   
     constexpr int N_READS = 4;
     constexpr int BM = 32;
     constexpr int BN = 32;
     col_reduce_looped<float, 2, BM, BN, N_READS><<<2, BM * BN / N_READS>>>(a, args);
     checkCuda(cudaGetLastError(), "kernel launch");
     checkCuda(cudaDeviceSynchronize(), "cudaDeviceSynchronize");
   
     cudaFree(a);
     printf("OK\n");
     return 0;
   }
3. zcbenz

   closed this as completedin #3273on Mar 18, 2026