|
| 1 | +#include "argmax.cuh" |
| 2 | +#include "utils.cuh" |
| 3 | +#include <cstddef> |
| 4 | +#include <cuda_bf16.hpp> |
| 5 | +#include <limits> |
| 6 | + |
| 7 | +namespace tensor::kernels { |
| 8 | + |
| 9 | +using namespace dtype; |
| 10 | + |
| 11 | +const int blockThreads = 256; |
| 12 | + |
| 13 | +__global__ void argmax_bfloat16_kernel(Cuda<int>* out, Cuda<bfloat16>* input, size_t num_reductions, size_t reduce_size, size_t reduce_stride) { |
| 14 | + __shared__ Cuda<bfloat16> shmem_val[blockThreads]; // NOLINT |
| 15 | + __shared__ int shmem_idx[blockThreads]; // NOLINT |
| 16 | + |
| 17 | + size_t tid = threadIdx.x; |
| 18 | + |
| 19 | + size_t reduction_idx = blockIdx.x; // which reduction are we doing? |
| 20 | + |
| 21 | + // decompose into outer and inner indices |
| 22 | + size_t outer_idx = reduction_idx / reduce_stride; |
| 23 | + size_t inner_idx = reduction_idx % reduce_stride; |
| 24 | + |
| 25 | + // base pointer for this reduction |
| 26 | + size_t base = (outer_idx * reduce_size * reduce_stride) + inner_idx; |
| 27 | + |
| 28 | + // reduce with a grid stride loop to handle reduce_size > blockThreads |
| 29 | + Cuda<bfloat16> thread_max_val = -std::numeric_limits<float>::infinity(); |
| 30 | + |
| 31 | + int thread_max_idx = tid; // NOLINT |
| 32 | + |
| 33 | + for (int element = tid; element < reduce_size; element += blockDim.x) { // NOLINT |
| 34 | + Cuda<bfloat16> incoming_val = input[base + (element * reduce_stride)]; |
| 35 | + int incoming_idx = element; |
| 36 | + |
| 37 | + if (incoming_val > thread_max_val) { |
| 38 | + thread_max_idx = incoming_idx; |
| 39 | + } |
| 40 | + |
| 41 | + thread_max_val = __hmax(thread_max_val, incoming_val); |
| 42 | + } |
| 43 | + |
| 44 | + // now we only have to reduce 'blockThreads' elements, which is easy within a block |
| 45 | + |
| 46 | + // load partial argmaxs onto shmem |
| 47 | + shmem_val[tid] = thread_max_val; |
| 48 | + shmem_idx[tid] = thread_max_idx; |
| 49 | + __syncthreads(); |
| 50 | + |
| 51 | + // reduce in shared memory |
| 52 | + for (int stride = blockDim.x / 2; stride > 32; stride >>= 1) { // NOLINT |
| 53 | + if (tid < stride) { |
| 54 | + Cuda<bfloat16> existing_val = shmem_val[tid]; |
| 55 | + Cuda<bfloat16> incoming_val = shmem_val[tid + stride]; |
| 56 | + int incoming_idx = shmem_idx[tid + stride]; |
| 57 | + if (incoming_val > existing_val) { |
| 58 | + shmem_idx[tid] = incoming_idx; |
| 59 | + } |
| 60 | + |
| 61 | + shmem_val[tid] = __hmax(existing_val, incoming_val); |
| 62 | + } |
| 63 | + __syncthreads(); |
| 64 | + } |
| 65 | + |
| 66 | + // warp shuffle for the final warp-level reduction |
| 67 | + if (tid < 32) { |
| 68 | + Cuda<bfloat16> val = shmem_val[tid]; |
| 69 | + int idx = shmem_idx[tid]; |
| 70 | + |
| 71 | + Cuda<bfloat16> incoming_val = shmem_val[tid + 32]; |
| 72 | + int incoming_idx = shmem_idx[tid + 32]; |
| 73 | + if (incoming_val > val) { |
| 74 | + shmem_idx[tid] = incoming_idx; |
| 75 | + } |
| 76 | + val = __hmax(val, incoming_val); |
| 77 | + |
| 78 | + for (int offset = 16; offset > 0; offset >>= 1) { |
| 79 | + Cuda<bfloat16> incoming_val = __shfl_down_sync(0xffffffff, val, offset); |
| 80 | + int incoming_idx = __shfl_down_sync(0xffffffff, idx, offset); |
| 81 | + if (incoming_val > val) { |
| 82 | + idx = incoming_idx; |
| 83 | + } |
| 84 | + val = __hmax(val, incoming_val); |
| 85 | + } |
| 86 | + |
| 87 | + if (tid == 0) { |
| 88 | + out[reduction_idx] = idx; |
| 89 | + } |
| 90 | + } |
| 91 | +} |
| 92 | + |
| 93 | +Tensor<int, CUDA> argmax_bfloat16(const TensorView<bfloat16, CUDA>& input, int dim, bool keepdim) { |
| 94 | + assert(input.is_contiguous() && "the tensor should be contiguous"); |
| 95 | + |
| 96 | + auto shape = input.shape; |
| 97 | + |
| 98 | + if (dim < 0) { |
| 99 | + dim = shape.size() + dim; |
| 100 | + } |
| 101 | + |
| 102 | + assert(dim >= 0 && static_cast<size_t>(dim) < shape.size()); |
| 103 | + |
| 104 | + size_t outer_size = 1; // how many reductions will we perform? ("batch size") |
| 105 | + size_t inner_size = 1; // what's the distance between elements to reduce? |
| 106 | + size_t reduce_size = 1; // how many elements each reduction needs to reduce over |
| 107 | + |
| 108 | + bool found_dim = false; |
| 109 | + |
| 110 | + // Output shape |
| 111 | + Shape out_shape; |
| 112 | + for (size_t i = 0; i < shape.size(); ++i) { |
| 113 | + if (i == static_cast<size_t>(dim)) { |
| 114 | + if (keepdim) { |
| 115 | + out_shape.push_back(1); |
| 116 | + } |
| 117 | + reduce_size = shape[dim]; |
| 118 | + found_dim = true; |
| 119 | + } else { |
| 120 | + if (!found_dim) { |
| 121 | + outer_size *= shape[i]; |
| 122 | + } else { |
| 123 | + inner_size *= shape[i]; |
| 124 | + } |
| 125 | + |
| 126 | + out_shape.push_back(shape[i]); |
| 127 | + } |
| 128 | + } |
| 129 | + |
| 130 | + if (out_shape.empty()) { |
| 131 | + out_shape.push_back(1); |
| 132 | + } |
| 133 | + |
| 134 | + auto n_elements = outer_size * inner_size; |
| 135 | + |
| 136 | + auto input_strides = get_all_strides(shape); |
| 137 | + |
| 138 | + TensorStorage<int, CUDA> storage(n_elements); |
| 139 | + Tensor<int, CUDA> out{out_shape, std::move(storage)}; |
| 140 | + |
| 141 | + int block_size = blockThreads; |
| 142 | + |
| 143 | + // Convert to device-native types for kernel call |
| 144 | + auto* out_d = reinterpret_cast<int*>(out.data()); // NOLINT |
| 145 | + auto* input_d = reinterpret_cast<Cuda<bfloat16>*>(input.data); // NOLINT |
| 146 | + |
| 147 | + argmax_bfloat16_kernel<<<n_elements, block_size>>>(out_d, input_d, n_elements, reduce_size, inner_size); |
| 148 | + |
| 149 | + return out; |
| 150 | +} |
| 151 | + |
| 152 | +} // namespace tensor::kernels |
0 commit comments