| 1 | /* |
| 2 | * Persistent TopK Scheduler for DSA Indexer |
| 3 | */ |
| 4 | |
| 5 | #ifndef PERSISTENT_TOPK_CUH_ |
| 6 | #define PERSISTENT_TOPK_CUH_ |
| 7 | |
| 8 | #include <cuda.h> |
| 9 | #include <cuda_fp16.h> |
| 10 | #include <cuda_runtime.h> |
| 11 | #include <cub/cub.cuh> |
| 12 | #include <cstdint> |
| 13 | |
| 14 | #include "topk_histogram_4096.cuh" |
| 15 | |
| 16 | namespace vllm { |
| 17 | namespace persistent { |
| 18 | |
| 19 | // ============================================================================ |
| 20 | // Constants |
| 21 | // ============================================================================ |
| 22 | |
| 23 | constexpr int kThreadsPerBlock = 1024; |
| 24 | constexpr int RADIX = 256; |
| 25 | |
| 26 | // Medium path: all shared state in dynamic smem (no static __shared__, |
| 27 | // which would inflate the kernel's smem footprint and kill occupancy |
| 28 | // for the decode/trivial paths). |
| 29 | constexpr size_t kMediumHistBytes = 2 * (RADIX + 128) * sizeof(int); // 3072 |
| 30 | constexpr size_t kMediumScalarsBytes = 5 * sizeof(int); // 20 |
| 31 | constexpr size_t kMediumHeaderSize = |
| 32 | (kMediumHistBytes + kMediumScalarsBytes + 127) & ~size_t(127); // 3200 |
| 33 | constexpr int MAX_BUFFERED_ITEMS = 4096; |
| 34 | constexpr size_t kSmemMedium = |
| 35 | kMediumHeaderSize + 2 * MAX_BUFFERED_ITEMS * sizeof(int); // 35968 |
| 36 | constexpr uint32_t RADIX_THRESHOLD = 32768; |
| 37 | |
| 38 | // Decode path constants |
| 39 | constexpr int kDecodeBins = 2048; |
| 40 | constexpr uint32_t HIST2048_THRESHOLD = 8192; |
| 41 | |
| 42 | // Large path: fixed shared memory for histograms + scalars |
| 43 | constexpr size_t kFixedSmemLarge = |
| 44 | ((RADIX + RADIX + 5) * sizeof(uint32_t) + 15) & ~size_t(15); |
| 45 | |
| 46 | // ============================================================================ |
| 47 | // Common helpers |
| 48 | // ============================================================================ |
| 49 | |
| 50 | __device__ __forceinline__ auto convert_to_uint32_v2(float x) -> uint32_t { |
| 51 | uint32_t bits = __float_as_uint(x); |
| 52 | return (bits & 0x80000000u) ? ~bits : (bits | 0x80000000u); |
| 53 | } |
| 54 | |
| 55 | __device__ __forceinline__ auto convert_to_uint8(float x) -> uint8_t { |
| 56 | __half h = __float2half_rn(x); |
| 57 | uint16_t bits = __half_as_ushort(h); |
| 58 | uint16_t key = (bits & 0x8000) ? static_cast<uint16_t>(~bits) |
| 59 | : static_cast<uint16_t>(bits | 0x8000); |
| 60 | return static_cast<uint8_t>(key >> 8); |
| 61 | } |
| 62 | |
| 63 | // ============================================================================ |
| 64 | // Vectorized load helpers |
| 65 | // ============================================================================ |
| 66 | |
| 67 | // Unconditional float4 load with cache hint (.cg = cache at global level only). |
| 68 | __device__ __forceinline__ void load_float4(const float* ptr, float& v0, |
| 69 | float& v1, float& v2, float& v3) { |
| 70 | uint32_t r0, r1, r2, r3; |
| 71 | asm volatile("ld.global.cg.v4.u32 {%0,%1,%2,%3}, [%4];\n" |
| 72 | : "=r"(r0), "=r"(r1), "=r"(r2), "=r"(r3) |
| 73 | : "l"(ptr)); |
| 74 | v0 = __uint_as_float(r0); |
| 75 | v1 = __uint_as_float(r1); |
| 76 | v2 = __uint_as_float(r2); |
| 77 | v3 = __uint_as_float(r3); |
| 78 | } |
| 79 | |
| 80 | // Per-element predicated scalar loads with -inf default. |
| 81 | __device__ __forceinline__ void load_float4_predicated(const float* ptr, |
| 82 | int base, int seq_len, |
| 83 | float& v0, float& v1, |
| 84 | float& v2, float& v3) { |
| 85 | uint32_t r0, r1, r2, r3; |
| 86 | int p0 = (base < seq_len); |
| 87 | int p1 = (base + 1 < seq_len); |
| 88 | int p2 = (base + 2 < seq_len); |
| 89 | int p3 = (base + 3 < seq_len); |
| 90 | asm volatile( |
| 91 | "{\n" |
| 92 | " .reg .pred pr0, pr1, pr2, pr3;\n" |
| 93 | " setp.ne.u32 pr0, %4, 0;\n" |
| 94 | " setp.ne.u32 pr1, %5, 0;\n" |
| 95 | " setp.ne.u32 pr2, %6, 0;\n" |
| 96 | " setp.ne.u32 pr3, %7, 0;\n" |
| 97 | " mov.u32 %0, 0xFF800000;\n" |
| 98 | " mov.u32 %1, 0xFF800000;\n" |
| 99 | " mov.u32 %2, 0xFF800000;\n" |
| 100 | " mov.u32 %3, 0xFF800000;\n" |
| 101 | " @pr0 ld.global.cg.u32 %0, [%8];\n" |
| 102 | " @pr1 ld.global.cg.u32 %1, [%8+4];\n" |
| 103 | " @pr2 ld.global.cg.u32 %2, [%8+8];\n" |
| 104 | " @pr3 ld.global.cg.u32 %3, [%8+12];\n" |
| 105 | "}\n" |
| 106 | : "=r"(r0), "=r"(r1), "=r"(r2), "=r"(r3) |
| 107 | : "r"(p0), "r"(p1), "r"(p2), "r"(p3), "l"(ptr)); |
| 108 | v0 = __uint_as_float(r0); |
| 109 | v1 = __uint_as_float(r1); |
| 110 | v2 = __uint_as_float(r2); |
| 111 | v3 = __uint_as_float(r3); |
| 112 | } |
| 113 | |
| 114 | // ============================================================================ |
| 115 | // Large path: inter-CTA coordination state (one per group) |
| 116 | // ============================================================================ |
| 117 | |
| 118 | struct RadixRowState { |
| 119 | uint32_t histogram[3][256]; // Triple-buffered histograms |
| 120 | uint32_t remaining_k; |
| 121 | uint32_t prefix; |
| 122 | int arrival_counter; |
| 123 | int output_counter; |
| 124 | }; |
| 125 | |
| 126 | // ============================================================================ |
| 127 | // Kernel parameters |
| 128 | // ============================================================================ |
| 129 | |
| 130 | struct PersistentTopKParams { |
| 131 | const float* __restrict__ input; // [num_rows, stride] |
| 132 | int32_t* __restrict__ output; // [num_rows, top_k] |
| 133 | const int32_t* __restrict__ lengths; // [num_rows] |
| 134 | RadixRowState* row_states; // large path: per-group state |
| 135 | uint32_t num_rows; |
| 136 | uint32_t stride; |
| 137 | uint32_t top_k; // actual k value for output stride |
| 138 | uint32_t chunk_size; // large path: elements per CTA |
| 139 | uint32_t ctas_per_group; // 1=medium, >1=large |
| 140 | uint32_t max_seq_len; // max seq_len across all rows (for early CTA exit) |
| 141 | }; |
| 142 | |
| 143 | // ============================================================================ |
| 144 | // Decode path: 2048-bin histogram for short sequences (seq_len <= 8192) |
| 145 | // Uses 11-bit half-precision bins for fine granularity. |
| 146 | // One histogram pass typically suffices since 8192/2048 = 4 elements/bin avg. |
| 147 | // ============================================================================ |
| 148 | |
| 149 | // 11-bit bin from half-precision representation (ascending: high values -> high |
| 150 | // bins) |
| 151 | __device__ __forceinline__ uint32_t decode_bin(float x) { |
| 152 | __half hx = __float2half(x); |
| 153 | uint16_t bits = __half_as_ushort(hx); |
| 154 | uint16_t key = (bits & 0x8000) ? static_cast<uint16_t>(~bits) |
| 155 | : static_cast<uint16_t>(bits | 0x8000); |
| 156 | return key >> 5; |
| 157 | } |
| 158 | |
| 159 | template <int TopK> |
| 160 | __device__ __noinline__ void histogram_2048_topk( |
| 161 | const float* __restrict__ logits, int32_t* __restrict__ output_indices, |
| 162 | int32_t seq_len) { |
| 163 | extern __shared__ int decode_smem[]; |
| 164 | const int tx = threadIdx.x; |
| 165 | const int lane = tx & 31; |
| 166 | |
| 167 | // ---- Layout constants ---- |
| 168 | constexpr int SBASE = 8192 - 8; // 8184 |
| 169 | constexpr int RHIST = RADIX + 128; // 384 |
| 170 | constexpr int BOFF = 2 * RHIST; // 768 |
| 171 | constexpr int DBUF = (SBASE - BOFF) / 2; // 3708 |
| 172 | constexpr int MAX_ITEMS_PER_THREAD = |
| 173 | (HIST2048_THRESHOLD + kThreadsPerBlock - 1) / kThreadsPerBlock; |
| 174 | |
| 175 | enum : int { sTHR = 0, sOUT = 1, sREF = 2, sFIN = 3, sBUF0 = 4, sBUF1 = 5 }; |
| 176 | |
| 177 | // ---- Initialize scalars (prevents stale data from prior rows) ---- |
| 178 | if (tx < 8) { |
| 179 | decode_smem[SBASE + tx] = 0; |
| 180 | } |
| 181 | |
| 182 | // ---- Phase 1: Build 2048-bin histogram with float4 vectorized loads ---- |
| 183 | int* histo = decode_smem; |
| 184 | uint16_t reg_bins[MAX_ITEMS_PER_THREAD]; |
| 185 | int nitems = 0; |
| 186 | |
| 187 | for (int i = tx; i < kDecodeBins; i += kThreadsPerBlock) { |
| 188 | histo[i] = 0; |
| 189 | } |
| 190 | __syncthreads(); |
| 191 | |
| 192 | const int n_vec = (seq_len + 3) >> 2; |
| 193 | const bool row_aligned = ((reinterpret_cast<uintptr_t>(logits) & 15) == 0); |
| 194 | |
| 195 | for (int i = tx; i < n_vec; i += kThreadsPerBlock) { |
| 196 | const int base = i << 2; |
| 197 | float v0, v1, v2, v3; |
| 198 | |
| 199 | if (row_aligned && base + 3 < seq_len) { |
| 200 | load_float4(logits + base, v0, v1, v2, v3); |
| 201 | } else { |
| 202 | load_float4_predicated(logits + base, base, seq_len, v0, v1, v2, v3); |
| 203 | } |
| 204 | |
| 205 | const uint16_t b0 = static_cast<uint16_t>(decode_bin(v0)); |
| 206 | const uint16_t b1 = static_cast<uint16_t>(decode_bin(v1)); |
| 207 | const uint16_t b2 = static_cast<uint16_t>(decode_bin(v2)); |
| 208 | const uint16_t b3 = static_cast<uint16_t>(decode_bin(v3)); |
| 209 | reg_bins[nitems++] = b0; |
| 210 | reg_bins[nitems++] = b1; |
| 211 | reg_bins[nitems++] = b2; |
| 212 | reg_bins[nitems++] = b3; |
| 213 | atomicAdd(&histo[b0], 1); |
| 214 | atomicAdd(&histo[b1], 1); |
| 215 | atomicAdd(&histo[b2], 1); |
| 216 | atomicAdd(&histo[b3], 1); |
| 217 | } |
| 218 | __syncthreads(); |
| 219 | |
| 220 | // ---- CUB suffix sum ---- |
| 221 | using BlockScanT = cub::BlockScan<int, kThreadsPerBlock>; |
| 222 | const int h0 = histo[2 * tx]; |
| 223 | const int pair_sum = h0 + histo[2 * tx + 1]; |
| 224 | |
| 225 | auto& scan_storage = *reinterpret_cast<typename BlockScanT::TempStorage*>( |
| 226 | decode_smem + kDecodeBins); |
| 227 | |
| 228 | int pair_prefix, total; |
| 229 | BlockScanT(scan_storage).ExclusiveSum(pair_sum, pair_prefix, total); |
| 230 | |
| 231 | // Find threshold bin purely from registers |
| 232 | const int pair_suffix = total - pair_prefix; |
| 233 | |
| 234 | if (pair_suffix >= TopK && (pair_suffix - h0) < TopK) { |
| 235 | decode_smem[SBASE + sTHR] = 2 * tx; |
| 236 | } |
| 237 | { |
| 238 | const int right_suf = pair_suffix - h0; |
| 239 | const int next_suf = pair_suffix - pair_sum; |
| 240 | if (right_suf >= TopK && next_suf < TopK) { |
| 241 | decode_smem[SBASE + sTHR] = 2 * tx + 1; |
| 242 | } |
| 243 | } |
| 244 | __syncthreads(); |
| 245 | |
| 246 | const int threshold = decode_smem[SBASE + sTHR]; |
| 247 | |
| 248 | // ---- Phase 2: Collection with warp-aggregated atomicAdds ---- |
| 249 | int* bufs[2] = {decode_smem + BOFF, decode_smem + BOFF + DBUF}; |
| 250 | const int sOUT_abs = SBASE + sOUT; |
| 251 | const int sBUF0_abs = SBASE + sBUF0; |
| 252 | |
| 253 | { |
| 254 | const uint32_t uthr = static_cast<uint32_t>(threshold); |
| 255 | int item = 0; |
| 256 | const int n_vec_iters = (n_vec + kThreadsPerBlock - 1) / kThreadsPerBlock; |
| 257 | |
| 258 | for (int iter = 0; iter < n_vec_iters; iter++) { |
| 259 | const int i = tx + iter * kThreadsPerBlock; |
| 260 | const bool vec_valid = (i < n_vec); |
| 261 | const int base_idx = i << 2; |
| 262 | |
| 263 | #pragma unroll 4 |
| 264 | for (int sub = 0; sub < 4; sub++) { |
| 265 | const int elem_idx = base_idx + sub; |
| 266 | uint32_t bin = 0; |
| 267 | if (vec_valid) bin = reg_bins[item++]; |
| 268 | const bool is_above = vec_valid && (bin > uthr); |
| 269 | const bool is_equal = vec_valid && (bin == uthr); |
| 270 | |
| 271 | const uint32_t above_mask = __ballot_sync(0xffffffff, is_above); |
| 272 | if (above_mask) { |
| 273 | const int above_count = __popc(above_mask); |
| 274 | const int above_rank = __popc(above_mask & ((1u << lane) - 1)); |
| 275 | int above_base; |
| 276 | if (lane == 0) { |
| 277 | above_base = atomicAdd(&decode_smem[sOUT_abs], above_count); |
| 278 | } |
| 279 | above_base = __shfl_sync(0xffffffff, above_base, 0); |
| 280 | if (is_above) { |
| 281 | output_indices[above_base + above_rank] = elem_idx; |
| 282 | } |
| 283 | } |
| 284 | |
| 285 | const uint32_t equal_mask = __ballot_sync(0xffffffff, is_equal); |
| 286 | if (equal_mask) { |
| 287 | const int equal_count = __popc(equal_mask); |
| 288 | const int equal_rank = __popc(equal_mask & ((1u << lane) - 1)); |
| 289 | int equal_base; |
| 290 | if (lane == 0) { |
| 291 | equal_base = atomicAdd(&decode_smem[sBUF0_abs], equal_count); |
| 292 | } |
| 293 | equal_base = __shfl_sync(0xffffffff, equal_base, 0); |
| 294 | if (is_equal && __builtin_expect(equal_base + equal_rank < DBUF, 1)) { |
| 295 | bufs[0][equal_base + equal_rank] = elem_idx; |
| 296 | } |
| 297 | } |
| 298 | } |
| 299 | } |
| 300 | } |
| 301 | __syncthreads(); |
| 302 | |
| 303 | int remaining_k = TopK - decode_smem[SBASE + sOUT]; |
| 304 | if (remaining_k <= 0) return; |
| 305 | |
| 306 | // If all buffered elements fit, output them all (common for short seqs) |
| 307 | const int raw_buf0 = decode_smem[SBASE + sBUF0]; |
| 308 | if (raw_buf0 <= remaining_k) { |
| 309 | const int nb = (raw_buf0 < DBUF) ? raw_buf0 : DBUF; |
| 310 | const int base = decode_smem[SBASE + sOUT]; |
| 311 | for (int i = tx; i < nb; i += kThreadsPerBlock) { |
| 312 | output_indices[base + i] = bufs[0][i]; |
| 313 | } |
| 314 | __syncthreads(); |
| 315 | return; |
| 316 | } |
| 317 | |
| 318 | // ---- Phase 3: Deferred refinement (rare path) ---- |
| 319 | int* refine[2] = {decode_smem, decode_smem + RHIST}; |
| 320 | const int num_buf0 = (raw_buf0 < DBUF) ? raw_buf0 : DBUF; |
| 321 | |
| 322 | for (int i = tx; i < RHIST; i += kThreadsPerBlock) { |
| 323 | refine[0][i] = 0; |
| 324 | } |
| 325 | __syncthreads(); |
| 326 | |
| 327 | for (int i = tx; i < num_buf0; i += kThreadsPerBlock) { |
| 328 | const uint32_t fp32 = convert_to_uint32_v2(logits[bufs[0][i]]); |
| 329 | atomicAdd(&refine[0][(fp32 >> 24) & 0xFF], 1); |
| 330 | } |
| 331 | __syncthreads(); |
| 332 | |
| 333 | auto compute_suffix_sum = [&]() { |
| 334 | #pragma unroll 8 |
| 335 | for (int i = 0; i < 8; ++i) { |
| 336 | if (tx < RADIX) { |
| 337 | const int stride = 1 << i; |
| 338 | const int s = i & 1; |
| 339 | const int d = s ^ 1; |
| 340 | int value = refine[s][tx]; |
| 341 | if (tx < RADIX - stride) value += refine[s][tx + stride]; |
| 342 | refine[d][tx] = value; |
| 343 | } |
| 344 | __syncthreads(); |
| 345 | } |
| 346 | }; |
| 347 | |
| 348 | #pragma unroll 4 |
| 349 | for (int pass = 0; pass < 4; ++pass) { |
| 350 | const int src = pass & 1; |
| 351 | const int dst = src ^ 1; |
| 352 | |
| 353 | const int raw_buf = decode_smem[SBASE + sBUF0 + src]; |
| 354 | const int num_buffered = (raw_buf < DBUF) ? raw_buf : DBUF; |
| 355 | |
| 356 | compute_suffix_sum(); |
| 357 | |
| 358 | if (tx < RADIX && refine[0][tx] > remaining_k && |
| 359 | refine[0][tx + 1] <= remaining_k) { |
| 360 | decode_smem[SBASE + sREF] = tx; |
| 361 | decode_smem[SBASE + sBUF0 + dst] = 0; |
| 362 | decode_smem[SBASE + sFIN] = remaining_k - refine[0][tx + 1]; |
| 363 | } |
| 364 | __syncthreads(); |
| 365 | |
| 366 | const int ref_thr = decode_smem[SBASE + sREF]; |
| 367 | remaining_k -= refine[0][ref_thr + 1]; |
| 368 | const int bit_offset = 24 - pass * 8; |
| 369 | |
| 370 | if (remaining_k == 0) { |
| 371 | for (int i = tx; i < num_buffered; i += kThreadsPerBlock) { |
| 372 | const int idx = bufs[src][i]; |
| 373 | const uint32_t fp32 = convert_to_uint32_v2(logits[idx]); |
| 374 | if (((fp32 >> bit_offset) & 0xFF) > static_cast<uint32_t>(ref_thr)) { |
| 375 | const int pos = atomicAdd(&decode_smem[SBASE + sOUT], 1); |
| 376 | output_indices[pos] = idx; |
| 377 | } |
| 378 | } |
| 379 | __syncthreads(); |
| 380 | break; |
| 381 | } |
| 382 | |
| 383 | __syncthreads(); |
| 384 | if (tx < RADIX + 1) refine[0][tx] = 0; |
| 385 | __syncthreads(); |
| 386 | |
| 387 | for (int i = tx; i < num_buffered; i += kThreadsPerBlock) { |
| 388 | const int idx = bufs[src][i]; |
| 389 | const float logit_val = logits[idx]; |
| 390 | const uint32_t fp32 = convert_to_uint32_v2(logit_val); |
| 391 | const int bin = (fp32 >> bit_offset) & 0xFF; |
| 392 | |
| 393 | if (bin > ref_thr) { |
| 394 | const int pos = atomicAdd(&decode_smem[SBASE + sOUT], 1); |
| 395 | output_indices[pos] = idx; |
| 396 | } else if (bin == ref_thr) { |
| 397 | if (pass == 3) { |
| 398 | const int slot = atomicAdd(&decode_smem[SBASE + sFIN], -1); |
| 399 | if (slot > 0) output_indices[TopK - slot] = idx; |
| 400 | } else { |
| 401 | const int bp = atomicAdd(&decode_smem[SBASE + sBUF0 + dst], 1); |
| 402 | if (__builtin_expect(bp < DBUF, 1)) { |
| 403 | bufs[dst][bp] = idx; |
| 404 | const int nbo = bit_offset - 8; |
| 405 | atomicAdd(&refine[0][(fp32 >> nbo) & 0xFF], 1); |
| 406 | } |
| 407 | } |
| 408 | } |
| 409 | } |
| 410 | __syncthreads(); |
| 411 | } |
| 412 | } |
| 413 | |
| 414 | // ============================================================================ |
| 415 | // Medium path: coarse FP16 histogram + 4-pass FP32 radix refinement |
| 416 | // For sequences 8K < seq_len <= 64K. |
| 417 | // ============================================================================ |
| 418 | |
| 419 | // Adapted from: |
| 420 | // https://github.com/sgl-project/sglang/blob/v0.5.8/sgl-kernel/csrc/elementwise/topk.cu#L87 |
| 421 | // by: DarkSharpness |
| 422 | // which at the same time is an optimized topk kernel copied from tilelang |
| 423 | // kernel |
| 424 | template <int TopK> |
| 425 | __device__ __noinline__ void histogram_256_topk( |
| 426 | const float* __restrict__ logits, int* __restrict__ output_indices, |
| 427 | int logits_offset, int seq_len) { |
| 428 | // All shared state lives in dynamic shared memory to avoid static |
| 429 | extern __shared__ char medium_smem[]; |
| 430 | |
| 431 | int (*shared_histogram)[RADIX + 128] = |
| 432 | reinterpret_cast<int (*)[RADIX + 128]>(medium_smem); |
| 433 | int* medium_scalars = reinterpret_cast<int*>(medium_smem + kMediumHistBytes); |
| 434 | int& shared_output_count = medium_scalars[0]; |
| 435 | int& shared_threshold_bin = medium_scalars[1]; |
| 436 | int* shared_buffered_count = &medium_scalars[2]; |
| 437 | int& shared_final_k = medium_scalars[4]; |
| 438 | int (*buffered_indices)[MAX_BUFFERED_ITEMS] = |
| 439 | reinterpret_cast<int (*)[MAX_BUFFERED_ITEMS]>(medium_smem + |
| 440 | kMediumHeaderSize); |
| 441 | |
| 442 | const int thread_id = threadIdx.x; |
| 443 | int remaining_k = TopK; |
| 444 | |
| 445 | if (thread_id < RADIX + 1) { |
| 446 | shared_histogram[0][thread_id] = 0; |
| 447 | } |
| 448 | __syncthreads(); |
| 449 | |
| 450 | for (int idx = thread_id; idx < seq_len; idx += kThreadsPerBlock) { |
| 451 | const auto bin = convert_to_uint8(logits[idx + logits_offset]); |
| 452 | atomicAdd(&shared_histogram[0][bin], 1); |
| 453 | } |
| 454 | __syncthreads(); |
| 455 | |
| 456 | auto compute_cumulative_sum = [&]() { |
| 457 | #pragma unroll 8 |
| 458 | for (int i = 0; i < 8; ++i) { |
| 459 | if (__builtin_expect(thread_id < RADIX, 1)) { |
| 460 | const int stride = 1 << i; |
| 461 | const int src_buffer = i & 1; |
| 462 | const int dst_buffer = src_buffer ^ 1; |
| 463 | int value = shared_histogram[src_buffer][thread_id]; |
| 464 | if (thread_id < RADIX - stride) { |
| 465 | value += shared_histogram[src_buffer][thread_id + stride]; |
| 466 | } |
| 467 | shared_histogram[dst_buffer][thread_id] = value; |
| 468 | } |
| 469 | __syncthreads(); |
| 470 | } |
| 471 | }; |
| 472 | |
| 473 | compute_cumulative_sum(); |
| 474 | |
| 475 | if (thread_id < RADIX && shared_histogram[0][thread_id] > remaining_k && |
| 476 | shared_histogram[0][thread_id + 1] <= remaining_k) { |
| 477 | shared_threshold_bin = thread_id; |
| 478 | shared_buffered_count[0] = 0; |
| 479 | shared_output_count = 0; |
| 480 | } |
| 481 | __syncthreads(); |
| 482 | |
| 483 | const int threshold_bin = shared_threshold_bin; |
| 484 | remaining_k -= shared_histogram[0][threshold_bin + 1]; |
| 485 | |
| 486 | if (remaining_k == 0) { |
| 487 | for (int idx = thread_id; idx < seq_len; idx += kThreadsPerBlock) { |
| 488 | const int bin = convert_to_uint8(logits[idx + logits_offset]); |
| 489 | if (bin > threshold_bin) { |
| 490 | const int output_pos = atomicAdd(&shared_output_count, 1); |
| 491 | output_indices[output_pos] = idx; |
| 492 | } |
| 493 | } |
| 494 | __syncthreads(); |
| 495 | return; |
| 496 | } |
| 497 | |
| 498 | __syncthreads(); |
| 499 | if (thread_id < RADIX + 1) { |
| 500 | shared_histogram[0][thread_id] = 0; |
| 501 | } |
| 502 | __syncthreads(); |
| 503 | |
| 504 | for (int idx = thread_id; idx < seq_len; idx += kThreadsPerBlock) { |
| 505 | const float logit_value = logits[idx + logits_offset]; |
| 506 | const int bin = convert_to_uint8(logit_value); |
| 507 | if (bin > threshold_bin) { |
| 508 | const int output_pos = atomicAdd(&shared_output_count, 1); |
| 509 | output_indices[output_pos] = idx; |
| 510 | } else if (bin == threshold_bin) { |
| 511 | const int buffer_pos = atomicAdd(&shared_buffered_count[0], 1); |
| 512 | if (__builtin_expect(buffer_pos < MAX_BUFFERED_ITEMS, 1)) { |
| 513 | buffered_indices[0][buffer_pos] = idx; |
| 514 | const uint32_t fp32_bits = convert_to_uint32_v2(logit_value); |
| 515 | const int next_bin = (fp32_bits >> 24) & 0xFF; |
| 516 | atomicAdd(&shared_histogram[0][next_bin], 1); |
| 517 | } |
| 518 | } |
| 519 | } |
| 520 | __syncthreads(); |
| 521 | |
| 522 | #pragma unroll 4 |
| 523 | for (int pass = 0; pass < 4; ++pass) { |
| 524 | const int src_buffer = pass % 2; |
| 525 | const int dst_buffer = src_buffer ^ 1; |
| 526 | const int raw_buffered = shared_buffered_count[src_buffer]; |
| 527 | const int num_buffered = |
| 528 | (raw_buffered < MAX_BUFFERED_ITEMS) ? raw_buffered : MAX_BUFFERED_ITEMS; |
| 529 | |
| 530 | compute_cumulative_sum(); |
| 531 | |
| 532 | if (thread_id < RADIX && shared_histogram[0][thread_id] > remaining_k && |
| 533 | shared_histogram[0][thread_id + 1] <= remaining_k) { |
| 534 | shared_threshold_bin = thread_id; |
| 535 | shared_buffered_count[dst_buffer] = 0; |
| 536 | shared_final_k = remaining_k - shared_histogram[0][thread_id + 1]; |
| 537 | } |
| 538 | __syncthreads(); |
| 539 | |
| 540 | const int threshold_bin = shared_threshold_bin; |
| 541 | remaining_k -= shared_histogram[0][threshold_bin + 1]; |
| 542 | const int bit_offset = 24 - pass * 8; |
| 543 | |
| 544 | if (remaining_k == 0) { |
| 545 | for (int i = thread_id; i < num_buffered; i += kThreadsPerBlock) { |
| 546 | const int idx = buffered_indices[src_buffer][i]; |
| 547 | const uint32_t fp32_bits = |
| 548 | convert_to_uint32_v2(logits[idx + logits_offset]); |
| 549 | const int bin = (fp32_bits >> bit_offset) & 0xFF; |
| 550 | if (bin > threshold_bin) { |
| 551 | const int output_pos = atomicAdd(&shared_output_count, 1); |
| 552 | output_indices[output_pos] = idx; |
| 553 | } |
| 554 | } |
| 555 | __syncthreads(); |
| 556 | break; |
| 557 | } |
| 558 | |
| 559 | __syncthreads(); |
| 560 | if (thread_id < RADIX + 1) { |
| 561 | shared_histogram[0][thread_id] = 0; |
| 562 | } |
| 563 | __syncthreads(); |
| 564 | |
| 565 | for (int i = thread_id; i < num_buffered; i += kThreadsPerBlock) { |
| 566 | const int idx = buffered_indices[src_buffer][i]; |
| 567 | const float logit_value = logits[idx + logits_offset]; |
| 568 | const uint32_t fp32_bits = convert_to_uint32_v2(logit_value); |
| 569 | const int bin = (fp32_bits >> bit_offset) & 0xFF; |
| 570 | if (bin > threshold_bin) { |
| 571 | const int output_pos = atomicAdd(&shared_output_count, 1); |
| 572 | output_indices[output_pos] = idx; |
| 573 | } else if (bin == threshold_bin) { |
| 574 | if (pass == 3) { |
| 575 | const int slot = atomicAdd(&shared_final_k, -1); |
| 576 | if (slot > 0) { |
| 577 | output_indices[TopK - slot] = idx; |
| 578 | } |
| 579 | } else { |
| 580 | const int buffer_pos = |
| 581 | atomicAdd(&shared_buffered_count[dst_buffer], 1); |
| 582 | if (__builtin_expect(buffer_pos < MAX_BUFFERED_ITEMS, 1)) { |
| 583 | buffered_indices[dst_buffer][buffer_pos] = idx; |
| 584 | const int next_bit_offset = bit_offset - 8; |
| 585 | const int next_bin = (fp32_bits >> next_bit_offset) & 0xFF; |
| 586 | atomicAdd(&shared_histogram[0][next_bin], 1); |
| 587 | } |
| 588 | } |
| 589 | } |
| 590 | } |
| 591 | __syncthreads(); |
| 592 | } |
| 593 | } |
| 594 | |
| 595 | // ============================================================================ |
| 596 | // Inter-CTA sync primitives |
| 597 | // ============================================================================ |
| 598 | |
| 599 | __device__ __forceinline__ int ld_acquire(int* ptr) { |
| 600 | int state = 0; |
| 601 | #if (__CUDA_ARCH__ >= 700) |
| 602 | asm volatile("ld.global.acquire.gpu.b32 %0, [%1];\n" |
| 603 | : "=r"(state) |
| 604 | : "l"(ptr)); |
| 605 | #else |
| 606 | asm volatile("ld.cg.global.b32 %0, [%1];\n" : "=r"(state) : "l"(ptr)); |
| 607 | #endif |
| 608 | return state; |
| 609 | } |
| 610 | |
| 611 | __device__ __forceinline__ void red_release(int* ptr, int val) { |
| 612 | #if (__CUDA_ARCH__ >= 700) |
| 613 | asm volatile("fence.acq_rel.gpu;\n"); |
| 614 | asm volatile("red.relaxed.gpu.global.add.s32 [%0], %1;\n" |
| 615 | : |
| 616 | : "l"(ptr), "r"(val)); |
| 617 | #else |
| 618 | __threadfence(); |
| 619 | atomicAdd(ptr, val); |
| 620 | #endif |
| 621 | } |
| 622 | |
| 623 | __device__ __forceinline__ void st_release(int* ptr, int val) { |
| 624 | #if (__CUDA_ARCH__ >= 700) |
| 625 | asm volatile("fence.acq_rel.gpu;\n"); |
| 626 | asm volatile("st.release.gpu.global.b32 [%0], %1;\n" : : "l"(ptr), "r"(val)); |
| 627 | #else |
| 628 | __threadfence(); |
| 629 | atomicExch(ptr, val); |
| 630 | #endif |
| 631 | } |
| 632 | |
| 633 | __device__ __forceinline__ void wait_ge(int* ptr, int target_val, |
| 634 | int thread_idx) { |
| 635 | if (thread_idx == 0) { |
| 636 | #pragma unroll 1 |
| 637 | while (ld_acquire(ptr) < target_val) { |
| 638 | } |
| 639 | } |
| 640 | __syncthreads(); |
| 641 | } |
| 642 | |
| 643 | // ============================================================================ |
| 644 | // Large path: multi-CTA radix select for sequences > 64K |
| 645 | // |
| 646 | // Each row is processed by a group of CTAs. Each CTA loads its chunk into |
| 647 | // shared memory as ordered uint32, then participates in 4 rounds of |
| 648 | // coordinated radix select via global-memory histograms and barriers. |
| 649 | // ============================================================================ |
| 650 | |
| 651 | // ============================================================================ |
| 652 | // Multi-CTA cooperative RadixTopK for a single large row. |
| 653 | // Adapted from https://github.com/flashinfer-ai/flashinfer/pull/2215 |
| 654 | // ============================================================================ |
| 655 | |
| 656 | template <int TopK, uint32_t VEC_SIZE> |
| 657 | __device__ void radix_topk(const float* __restrict__ row_input, |
| 658 | int32_t* __restrict__ row_output, uint32_t seq_len, |
| 659 | uint32_t my_chunk_start, uint32_t chunk_size, |
| 660 | uint32_t* local_histogram, uint32_t* suffix_sum, |
| 661 | uint32_t* shared_scalars, uint32_t* shared_ordered, |
| 662 | RadixRowState* state, uint32_t cta_in_group, |
| 663 | uint32_t ctas_per_group, int& barrier_phase, |
| 664 | uint32_t radix_iter, uint32_t tx) { |
| 665 | const uint32_t my_chunk_end = (my_chunk_start + chunk_size < seq_len) |
| 666 | ? my_chunk_start + chunk_size |
| 667 | : seq_len; |
| 668 | const uint32_t actual_chunk_size = |
| 669 | (my_chunk_start < seq_len) ? (my_chunk_end - my_chunk_start) : 0; |
| 670 | |
| 671 | // -- Stage 1: Load chunk to shared memory as ordered uint32 -- |
| 672 | { |
| 673 | const uint32_t aligned_size = (actual_chunk_size / VEC_SIZE) * VEC_SIZE; |
| 674 | |
| 675 | for (uint32_t i = tx * VEC_SIZE; i < aligned_size; |
| 676 | i += kThreadsPerBlock * VEC_SIZE) { |
| 677 | const float* src = row_input + my_chunk_start + i; |
| 678 | if constexpr (VEC_SIZE == 4) { |
| 679 | float4 v = *reinterpret_cast<const float4*>(src); |
| 680 | shared_ordered[i] = convert_to_uint32_v2(v.x); |
| 681 | shared_ordered[i + 1] = convert_to_uint32_v2(v.y); |
| 682 | shared_ordered[i + 2] = convert_to_uint32_v2(v.z); |
| 683 | shared_ordered[i + 3] = convert_to_uint32_v2(v.w); |
| 684 | } else if constexpr (VEC_SIZE == 2) { |
| 685 | float2 v = *reinterpret_cast<const float2*>(src); |
| 686 | shared_ordered[i] = convert_to_uint32_v2(v.x); |
| 687 | shared_ordered[i + 1] = convert_to_uint32_v2(v.y); |
| 688 | } else { |
| 689 | shared_ordered[i] = convert_to_uint32_v2(*src); |
| 690 | } |
| 691 | } |
| 692 | for (uint32_t i = aligned_size + tx; i < actual_chunk_size; |
| 693 | i += kThreadsPerBlock) { |
| 694 | shared_ordered[i] = convert_to_uint32_v2(row_input[my_chunk_start + i]); |
| 695 | } |
| 696 | } |
| 697 | __syncthreads(); |
| 698 | |
| 699 | // -- Init radix select state -- |
| 700 | if (tx == 0) { |
| 701 | shared_scalars[0] = 0; // prefix |
| 702 | shared_scalars[1] = TopK; // remaining_k |
| 703 | } |
| 704 | __syncthreads(); |
| 705 | |
| 706 | // -- Initial barrier -- |
| 707 | if (tx == 0) { |
| 708 | red_release(&state->arrival_counter, 1); |
| 709 | } |
| 710 | wait_ge(&state->arrival_counter, |
| 711 | (barrier_phase + 1) * static_cast<int>(ctas_per_group), tx); |
| 712 | barrier_phase++; |
| 713 | __syncthreads(); |
| 714 | |
| 715 | if (cta_in_group == 0 && tx == 0) { |
| 716 | st_release(&state->output_counter, 0); |
| 717 | } |
| 718 | |
| 719 | // -- Stage 2: 4 rounds of radix select -- |
| 720 | for (uint32_t round = 0; round < 4; round++) { |
| 721 | const uint32_t global_round = radix_iter * 4 + round; |
| 722 | const uint32_t shift = 24 - round * 8; |
| 723 | const uint32_t prefix = shared_scalars[0]; |
| 724 | const uint32_t remaining_k = shared_scalars[1]; |
| 725 | |
| 726 | uint32_t* current_hist = state->histogram[global_round % 3]; |
| 727 | uint32_t* next_hist = state->histogram[(global_round + 1) % 3]; |
| 728 | |
| 729 | for (uint32_t i = tx; i < RADIX; i += kThreadsPerBlock) { |
| 730 | local_histogram[i] = 0; |
| 731 | } |
| 732 | __syncthreads(); |
| 733 | |
| 734 | for (uint32_t i = tx; i < actual_chunk_size; i += kThreadsPerBlock) { |
| 735 | uint32_t ordered = shared_ordered[i]; |
| 736 | uint32_t mask = (round == 0) ? 0u : (~0u << (32 - round * 8)); |
| 737 | if ((ordered & mask) == prefix) { |
| 738 | uint32_t bucket = (ordered >> shift) & 0xFF; |
| 739 | atomicAdd(&local_histogram[bucket], 1); |
| 740 | } |
| 741 | } |
| 742 | __syncthreads(); |
| 743 | |
| 744 | for (uint32_t i = tx; i < RADIX; i += kThreadsPerBlock) { |
| 745 | if (local_histogram[i] > 0) { |
| 746 | atomicAdd(¤t_hist[i], local_histogram[i]); |
| 747 | } |
| 748 | } |
| 749 | |
| 750 | if (cta_in_group == 0) { |
| 751 | for (uint32_t i = tx; i < RADIX; i += kThreadsPerBlock) { |
| 752 | next_hist[i] = 0; |
| 753 | } |
| 754 | } |
| 755 | |
| 756 | if (tx == 0) { |
| 757 | red_release(&state->arrival_counter, 1); |
| 758 | } |
| 759 | wait_ge(&state->arrival_counter, |
| 760 | (barrier_phase + 1) * static_cast<int>(ctas_per_group), tx); |
| 761 | barrier_phase++; |
| 762 | __syncthreads(); |
| 763 | |
| 764 | for (uint32_t i = tx; i < RADIX; i += kThreadsPerBlock) { |
| 765 | suffix_sum[i] = current_hist[i]; |
| 766 | } |
| 767 | __syncthreads(); |
| 768 | |
| 769 | for (uint32_t stride = 1; stride < RADIX; stride *= 2) { |
| 770 | uint32_t val = 0; |
| 771 | if (tx < RADIX) { |
| 772 | val = suffix_sum[tx]; |
| 773 | if (tx + stride < RADIX) val += suffix_sum[tx + stride]; |
| 774 | } |
| 775 | __syncthreads(); |
| 776 | if (tx < RADIX) suffix_sum[tx] = val; |
| 777 | __syncthreads(); |
| 778 | } |
| 779 | |
| 780 | if (tx == 0) { |
| 781 | shared_scalars[2] = 0; |
| 782 | shared_scalars[3] = remaining_k; |
| 783 | } |
| 784 | __syncthreads(); |
| 785 | |
| 786 | if (tx < RADIX) { |
| 787 | uint32_t count_ge = suffix_sum[tx]; |
| 788 | uint32_t count_gt = (tx + 1 < RADIX) ? suffix_sum[tx + 1] : 0; |
| 789 | if (count_ge >= remaining_k && count_gt < remaining_k) { |
| 790 | shared_scalars[2] = tx; |
| 791 | shared_scalars[3] = remaining_k - count_gt; |
| 792 | } |
| 793 | } |
| 794 | __syncthreads(); |
| 795 | |
| 796 | if (tx == 0) { |
| 797 | shared_scalars[0] = prefix | (shared_scalars[2] << shift); |
| 798 | shared_scalars[1] = shared_scalars[3]; |
| 799 | } |
| 800 | __syncthreads(); |
| 801 | } // end 4 radix rounds |
| 802 | |
| 803 | // -- Count local > pivot elements -- |
| 804 | const uint32_t ordered_pivot = shared_scalars[0]; |
| 805 | |
| 806 | if (tx == 0) suffix_sum[0] = 0; |
| 807 | __syncthreads(); |
| 808 | |
| 809 | uint32_t my_gt_count = 0; |
| 810 | for (uint32_t i = tx; i < actual_chunk_size; i += kThreadsPerBlock) { |
| 811 | if (shared_ordered[i] > ordered_pivot) my_gt_count++; |
| 812 | } |
| 813 | for (int offset = 16; offset > 0; offset /= 2) { |
| 814 | my_gt_count += __shfl_down_sync(0xffffffff, my_gt_count, offset); |
| 815 | } |
| 816 | if (tx % 32 == 0 && my_gt_count > 0) { |
| 817 | atomicAdd(&suffix_sum[0], my_gt_count); |
| 818 | } |
| 819 | __syncthreads(); |
| 820 | const uint32_t local_gt_count = suffix_sum[0]; |
| 821 | |
| 822 | // -- Stage 3: Collect top-k indices -- |
| 823 | if (tx == 0) { |
| 824 | local_histogram[0] = 0; |
| 825 | if (local_gt_count > 0) { |
| 826 | local_histogram[1] = |
| 827 | atomicAdd(&state->output_counter, static_cast<int>(local_gt_count)); |
| 828 | } |
| 829 | } |
| 830 | __syncthreads(); |
| 831 | |
| 832 | for (uint32_t i = tx; i < actual_chunk_size; i += kThreadsPerBlock) { |
| 833 | if (shared_ordered[i] > ordered_pivot) { |
| 834 | uint32_t local_pos = atomicAdd(&local_histogram[0], 1); |
| 835 | int pos = static_cast<int>(local_histogram[1]) + local_pos; |
| 836 | row_output[pos] = static_cast<int32_t>(my_chunk_start + i); |
| 837 | } |
| 838 | } |
| 839 | |
| 840 | if (tx == 0) { |
| 841 | red_release(&state->arrival_counter, 1); |
| 842 | } |
| 843 | wait_ge(&state->arrival_counter, |
| 844 | (barrier_phase + 1) * static_cast<int>(ctas_per_group), tx); |
| 845 | barrier_phase++; |
| 846 | __syncthreads(); |
| 847 | |
| 848 | for (uint32_t i = tx; i < actual_chunk_size; i += kThreadsPerBlock) { |
| 849 | if (shared_ordered[i] == ordered_pivot) { |
| 850 | int pos = atomicAdd(&state->output_counter, 1); |
| 851 | if (pos < TopK) { |
| 852 | row_output[pos] = static_cast<int32_t>(my_chunk_start + i); |
| 853 | } |
| 854 | } |
| 855 | } |
| 856 | } |
| 857 | |
| 858 | // ============================================================================ |
| 859 | // Persistent kernel — BS≤32, decode/medium/large paths with RadixTopK |
| 860 | // BS>32 uses standalone histogram_256_buffered_topk (separate kernel, |
| 861 | // see filtered_topk.cuh) |
| 862 | // ============================================================================ |
| 863 | |
| 864 | template <int TopK = 2048, uint32_t VEC_SIZE = 1> |
| 865 | __global__ void __launch_bounds__(kThreadsPerBlock, 2) |
| 866 | persistent_topk_kernel(PersistentTopKParams params) { |
| 867 | const uint32_t tx = threadIdx.x; |
| 868 | extern __shared__ uint8_t smem_raw[]; |
| 869 | |
| 870 | // ======================================================================== |
| 871 | // Group mode: multi-CTA groups with static round-robin row assignment. |
| 872 | // Non-large rows: CTA-0 handles trivial/decode/medium. |
| 873 | // Large rows: all CTAs in the group cooperate via RadixTopK. |
| 874 | // ======================================================================== |
| 875 | const uint32_t ctas_per_group = params.ctas_per_group; |
| 876 | const uint32_t group_id = blockIdx.x / ctas_per_group; |
| 877 | const uint32_t cta_in_group = blockIdx.x % ctas_per_group; |
| 878 | const uint32_t num_groups = gridDim.x / ctas_per_group; |
| 879 | const uint32_t chunk_size = params.chunk_size; |
| 880 | |
| 881 | if (blockIdx.x >= num_groups * ctas_per_group) return; |
| 882 | |
| 883 | // Early exit: non-CTA-0 threads are never needed if no large rows exist |
| 884 | if (cta_in_group != 0 && params.max_seq_len <= RADIX_THRESHOLD) return; |
| 885 | |
| 886 | uint32_t* local_histogram = reinterpret_cast<uint32_t*>(smem_raw); |
| 887 | uint32_t* suffix_sum = local_histogram + RADIX; |
| 888 | uint32_t* shared_scalars = suffix_sum + RADIX; |
| 889 | uint32_t* shared_ordered = |
| 890 | reinterpret_cast<uint32_t*>(smem_raw + kFixedSmemLarge); |
| 891 | |
| 892 | // RadixRowState for multi-CTA cooperative radix. |
| 893 | // Zero-initialization is done host-side via cudaMemsetAsync in topk.cu |
| 894 | // before launch — that gives a stream-ordered happens-before edge for all |
| 895 | // CTAs, which the previous in-kernel init (CTA-0 only + intra-CTA |
| 896 | // __syncthreads) did not provide and which manifested as a race against |
| 897 | // CTA-1+'s first red_release on arrival_counter. |
| 898 | RadixRowState* state = ¶ms.row_states[group_id]; |
| 899 | |
| 900 | int barrier_phase = 0; |
| 901 | uint32_t radix_iter = 0; |
| 902 | const uint32_t total_iters = (params.num_rows + num_groups - 1) / num_groups; |
| 903 | |
| 904 | for (uint32_t iter = 0; iter < total_iters; iter++) { |
| 905 | // Static round-robin: all CTAs in the group implicitly agree on the row |
| 906 | uint32_t row_idx = group_id + iter * num_groups; |
| 907 | if (row_idx >= params.num_rows) break; |
| 908 | |
| 909 | // Clamp the row length before any decision is made on it. |
| 910 | // |
| 911 | // `lengths` is int32 and is consumed here as uint32, so a negative value |
| 912 | // (e.g. a padded decode slot whose per-token context length underflowed) |
| 913 | // would reinterpret as ~4e9 and sail past every threshold below. Any |
| 914 | // value beyond the row width would also read into the next row. |
| 915 | // |
| 916 | // Clamping to max_seq_len additionally keeps this per-row decision |
| 917 | // consistent with the `cta_in_group != 0` early exit above, which is |
| 918 | // taken from the host-side scalar: when max_seq_len <= RADIX_THRESHOLD |
| 919 | // the non-leader CTAs return immediately, so a leader that reached the |
| 920 | // cooperative radix path would wait on the inter-CTA barrier for peers |
| 921 | // that no longer exist and spin until the kernel is killed. |
| 922 | const int32_t raw_len = params.lengths[row_idx]; |
| 923 | const uint32_t row_bound = |
| 924 | params.stride < params.max_seq_len ? params.stride : params.max_seq_len; |
| 925 | const uint32_t non_negative_len = |
| 926 | raw_len > 0 ? static_cast<uint32_t>(raw_len) : 0u; |
| 927 | const uint32_t seq_len = |
| 928 | non_negative_len < row_bound ? non_negative_len : row_bound; |
| 929 | int32_t* row_output = params.output + row_idx * params.top_k; |
| 930 | const float* row_input = params.input + row_idx * params.stride; |
| 931 | |
| 932 | if (seq_len <= RADIX_THRESHOLD) { |
| 933 | if (cta_in_group == 0) { |
| 934 | if (seq_len <= static_cast<uint32_t>(TopK)) { |
| 935 | // Trivial case: seq_len <= TopK |
| 936 | for (uint32_t i = tx; i < static_cast<uint32_t>(TopK); |
| 937 | i += kThreadsPerBlock) { |
| 938 | row_output[i] = (i < seq_len) ? static_cast<int32_t>(i) : -1; |
| 939 | } |
| 940 | } else if (seq_len <= static_cast<uint32_t>(HIST2048_THRESHOLD)) { |
| 941 | histogram_2048_topk<TopK>(row_input, row_output, seq_len); |
| 942 | } else { |
| 943 | histogram_256_topk<TopK>(row_input, row_output, 0, seq_len); |
| 944 | } |
| 945 | } |
| 946 | continue; |
| 947 | } |
| 948 | |
| 949 | const uint32_t my_chunk_start = cta_in_group * chunk_size; |
| 950 | radix_topk<TopK, VEC_SIZE>( |
| 951 | row_input, row_output, seq_len, my_chunk_start, chunk_size, |
| 952 | local_histogram, suffix_sum, shared_scalars, shared_ordered, state, |
| 953 | cta_in_group, ctas_per_group, barrier_phase, radix_iter, tx); |
| 954 | radix_iter++; |
| 955 | } |
| 956 | } |
| 957 | |
| 958 | } // namespace persistent |
| 959 | |
| 960 | // ============================================================================ |
| 961 | // ============================================================================ |
| 962 | // Optimized FilteredTopK — single CTA per row for bs > 32. |
| 963 | // Kept with persistent_topk so the portable fallback owns the non-cluster path. |
| 964 | // ============================================================================ |
| 965 | namespace filtered_topk { |
| 966 | |
| 967 | namespace hist4096 = topk_histogram_4096; |
| 968 | |
| 969 | // ============================================================================ |
| 970 | // FilteredTopK — single CTA per row for bs > 32 |
| 971 | // Adapted from https://github.com/flashinfer-ai/flashinfer/pull/2215 |
| 972 | // ============================================================================ |
| 973 | |
| 974 | #define FLASHINFER_CUDA_CALL(func, ...) \ |
| 975 | { \ |
| 976 | cudaError_t e = (func); \ |
| 977 | if (e != cudaSuccess) { \ |
| 978 | return e; \ |
| 979 | } \ |
| 980 | } |
| 981 | |
| 982 | #define FLASHINFER_INLINE inline __attribute__((always_inline)) __device__ |
| 983 | |
| 984 | template <typename T, size_t N> |
| 985 | struct vec_t { |
| 986 | T data[N]; |
| 987 | |
| 988 | FLASHINFER_INLINE T& operator[](size_t i) { return data[i]; } |
| 989 | FLASHINFER_INLINE const T& operator[](size_t i) const { return data[i]; } |
| 990 | |
| 991 | FLASHINFER_INLINE void cast_load(const T* ptr) { |
| 992 | #pragma unroll |
| 993 | for (size_t i = 0; i < N; ++i) { |
| 994 | data[i] = ptr[i]; |
| 995 | } |
| 996 | } |
| 997 | }; |
| 998 | #undef FLASHINFER_INLINE |
| 999 | |
| 1000 | // FilteredTopK traits for different data types |
| 1001 | template <typename DType> |
| 1002 | struct FilteredTopKTraits; |
| 1003 | |
| 1004 | // Specialization for float (32-bit): coarse histogram uses FP16 high 8 bits, 4 |
| 1005 | // refinement rounds |
| 1006 | template <> |
| 1007 | struct FilteredTopKTraits<float> { |
| 1008 | using OrderedType = uint32_t; |
| 1009 | static constexpr int NUM_REFINE_ROUNDS = 4; |
| 1010 | static constexpr int FIRST_REFINE_SHIFT = 24; |
| 1011 | |
| 1012 | __device__ __forceinline__ static uint8_t ToCoarseKey(float x) { |
| 1013 | // Convert to FP16 representation and extract high 8 bits |
| 1014 | __half h = __float2half_rn(x); |
| 1015 | uint16_t bits = __half_as_ushort(h); |
| 1016 | uint16_t key = (bits & 0x8000) ? static_cast<uint16_t>(~bits) |
| 1017 | : static_cast<uint16_t>(bits | 0x8000); |
| 1018 | return static_cast<uint8_t>(key >> 8); |
| 1019 | } |
| 1020 | |
| 1021 | __device__ __forceinline__ static OrderedType ToOrdered(float x) { |
| 1022 | uint32_t bits = __float_as_uint(x); |
| 1023 | return (bits & 0x80000000u) ? ~bits : (bits | 0x80000000u); |
| 1024 | } |
| 1025 | }; |
| 1026 | |
| 1027 | constexpr uint32_t FILTERED_TOPK_BLOCK_THREADS = 1024; |
| 1028 | constexpr uint32_t FILTERED_TOPK_SMEM_INPUT_SIZE = |
| 1029 | 16 * 1024; // 16K indices per buffer |
| 1030 | constexpr size_t FILTERED_TOPK_SMEM_DYNAMIC = |
| 1031 | sizeof(int) * 2 * FILTERED_TOPK_SMEM_INPUT_SIZE; // 128KB |
| 1032 | |
| 1033 | template <uint32_t MAX_K> |
| 1034 | struct FilteredTopKStorage { |
| 1035 | alignas(128) int histogram[2][256 + 128]; |
| 1036 | alignas(128) int counter; |
| 1037 | alignas(128) int threshold_bin; |
| 1038 | alignas(128) int num_input[2]; |
| 1039 | alignas(128) int indices[MAX_K]; |
| 1040 | int last_remain; |
| 1041 | }; |
| 1042 | |
| 1043 | // With CheckOverflow, return false before using a truncated stash so the caller |
| 1044 | // can retry with exact full-row selection. |
| 1045 | template <typename DType, typename IdType, int VEC_SIZE, uint32_t MAX_K, |
| 1046 | bool UsePredicatedShortLoads, bool CheckOverflow = false> |
| 1047 | __device__ __forceinline__ bool filtered_topk_row( |
| 1048 | const DType* score, IdType* dst, int length, uint32_t top_k, |
| 1049 | FilteredTopKStorage<MAX_K>& storage) { |
| 1050 | constexpr uint32_t BLOCK_SIZE = FILTERED_TOPK_BLOCK_THREADS; |
| 1051 | constexpr int RADIX = 256; |
| 1052 | constexpr int SMEM_INPUT_SIZE = FILTERED_TOPK_SMEM_INPUT_SIZE; |
| 1053 | const int tx = threadIdx.x; |
| 1054 | |
| 1055 | // Trivial case: length <= top_k |
| 1056 | if (length <= static_cast<int>(top_k)) { |
| 1057 | for (int i = tx; i < static_cast<int>(top_k); i += BLOCK_SIZE) { |
| 1058 | dst[i] = (i < length) ? static_cast<IdType>(i) : static_cast<IdType>(-1); |
| 1059 | } |
| 1060 | return true; |
| 1061 | } |
| 1062 | |
| 1063 | // Short path |
| 1064 | if (length <= 32768) { |
| 1065 | extern __shared__ uint8_t _smem_reg[]; |
| 1066 | if constexpr (UsePredicatedShortLoads) { |
| 1067 | hist4096::histogram_4096_topk_predicated<MAX_K, 12, 8>(score, dst, length, |
| 1068 | _smem_reg); |
| 1069 | } else { |
| 1070 | hist4096::histogram_4096_topk<MAX_K, 12, 8>(score, dst, length, |
| 1071 | _smem_reg); |
| 1072 | } |
| 1073 | return true; |
| 1074 | } |
| 1075 | |
| 1076 | auto& s_histogram_buf = storage.histogram; |
| 1077 | auto& s_counter = storage.counter; |
| 1078 | auto& s_threshold_bin_id = storage.threshold_bin; |
| 1079 | auto& s_num_input = storage.num_input; |
| 1080 | auto& s_indices = storage.indices; |
| 1081 | auto& s_histogram = s_histogram_buf[0]; |
| 1082 | |
| 1083 | // Dynamic shared memory for input double buffer |
| 1084 | extern __shared__ int s_input_idx[][SMEM_INPUT_SIZE]; |
| 1085 | |
| 1086 | using Traits = FilteredTopKTraits<DType>; |
| 1087 | int topk = top_k; |
| 1088 | |
| 1089 | // Stage 1: 8-bit coarse histogram with vectorized loads |
| 1090 | if (tx < RADIX + 1) s_histogram[tx] = 0; |
| 1091 | __syncthreads(); |
| 1092 | |
| 1093 | vec_t<DType, VEC_SIZE> score_vec; |
| 1094 | |
| 1095 | const int aligned_length = (length / VEC_SIZE) * VEC_SIZE; |
| 1096 | #pragma unroll 2 |
| 1097 | for (int base = tx * VEC_SIZE; base < aligned_length; |
| 1098 | base += BLOCK_SIZE * VEC_SIZE) { |
| 1099 | score_vec.cast_load(&score[base]); |
| 1100 | #pragma unroll |
| 1101 | for (int j = 0; j < VEC_SIZE; ++j) { |
| 1102 | const auto bin = Traits::ToCoarseKey(score_vec[j]); |
| 1103 | atomicAdd(&s_histogram[bin], 1); |
| 1104 | } |
| 1105 | } |
| 1106 | // Handle tail |
| 1107 | for (int i = aligned_length + tx; i < length; i += BLOCK_SIZE) { |
| 1108 | const auto bin = Traits::ToCoarseKey(score[i]); |
| 1109 | atomicAdd(&s_histogram[bin], 1); |
| 1110 | } |
| 1111 | __syncthreads(); |
| 1112 | |
| 1113 | // Suffix sum |
| 1114 | const auto run_cumsum = [&]() { |
| 1115 | #pragma unroll 8 |
| 1116 | for (int i = 0; i < 8; ++i) { |
| 1117 | if (tx < RADIX) { |
| 1118 | const auto j = 1 << i; |
| 1119 | const auto k = i & 1; |
| 1120 | auto value = s_histogram_buf[k][tx]; |
| 1121 | if (tx < RADIX - j) { |
| 1122 | value += s_histogram_buf[k][tx + j]; |
| 1123 | } |
| 1124 | s_histogram_buf[k ^ 1][tx] = value; |
| 1125 | } |
| 1126 | __syncthreads(); |
| 1127 | } |
| 1128 | }; |
| 1129 | |
| 1130 | run_cumsum(); |
| 1131 | if (tx < RADIX && s_histogram[tx] > topk && s_histogram[tx + 1] <= topk) { |
| 1132 | s_threshold_bin_id = tx; |
| 1133 | s_num_input[0] = 0; |
| 1134 | s_counter = 0; |
| 1135 | } |
| 1136 | __syncthreads(); |
| 1137 | |
| 1138 | const auto threshold_bin = s_threshold_bin_id; |
| 1139 | topk -= s_histogram[threshold_bin + 1]; |
| 1140 | |
| 1141 | constexpr int NUM_ROUNDS = Traits::NUM_REFINE_ROUNDS; |
| 1142 | constexpr int FIRST_SHIFT = Traits::FIRST_REFINE_SHIFT; |
| 1143 | |
| 1144 | if (topk == 0) { |
| 1145 | // Collect indices where bin > threshold |
| 1146 | #pragma unroll 2 |
| 1147 | for (int base = tx * VEC_SIZE; base < aligned_length; |
| 1148 | base += BLOCK_SIZE * VEC_SIZE) { |
| 1149 | score_vec.cast_load(&score[base]); |
| 1150 | #pragma unroll |
| 1151 | for (int j = 0; j < VEC_SIZE; ++j) { |
| 1152 | const auto bin = static_cast<int>(Traits::ToCoarseKey(score_vec[j])); |
| 1153 | if (bin > threshold_bin) { |
| 1154 | const auto pos = atomicAdd(&s_counter, 1); |
| 1155 | s_indices[pos] = base + j; |
| 1156 | } |
| 1157 | } |
| 1158 | } |
| 1159 | // Handle tail |
| 1160 | for (int i = aligned_length + tx; i < length; i += BLOCK_SIZE) { |
| 1161 | const auto bin = static_cast<int>(Traits::ToCoarseKey(score[i])); |
| 1162 | if (bin > threshold_bin) { |
| 1163 | const auto pos = atomicAdd(&s_counter, 1); |
| 1164 | s_indices[pos] = i; |
| 1165 | } |
| 1166 | } |
| 1167 | __syncthreads(); |
| 1168 | } else { |
| 1169 | __syncthreads(); |
| 1170 | if (tx < RADIX + 1) s_histogram[tx] = 0; |
| 1171 | __syncthreads(); |
| 1172 | |
| 1173 | // Filter + histogram for refinement |
| 1174 | auto filter_and_add_to_histogram = [&](auto raw_input, int index) { |
| 1175 | const auto bin = static_cast<int>(Traits::ToCoarseKey(raw_input)); |
| 1176 | if (bin > threshold_bin) { |
| 1177 | const auto pos = atomicAdd(&s_counter, 1); |
| 1178 | s_indices[pos] = index; |
| 1179 | } else if (bin == threshold_bin) { |
| 1180 | const auto pos = atomicAdd(&s_num_input[0], 1); |
| 1181 | if (__builtin_expect(pos < SMEM_INPUT_SIZE, 1)) { |
| 1182 | s_input_idx[0][pos] = index; |
| 1183 | const auto ordered = Traits::ToOrdered(raw_input); |
| 1184 | const auto sub_bin = (ordered >> FIRST_SHIFT) & 0xFF; |
| 1185 | atomicAdd(&s_histogram[sub_bin], 1); |
| 1186 | } |
| 1187 | } |
| 1188 | }; |
| 1189 | #pragma unroll 2 |
| 1190 | for (int base = tx * VEC_SIZE; base < aligned_length; |
| 1191 | base += BLOCK_SIZE * VEC_SIZE) { |
| 1192 | score_vec.cast_load(&score[base]); |
| 1193 | #pragma unroll |
| 1194 | for (int j = 0; j < VEC_SIZE; ++j) { |
| 1195 | filter_and_add_to_histogram(score_vec[j], base + j); |
| 1196 | } |
| 1197 | } |
| 1198 | // Handle tail |
| 1199 | for (int i = aligned_length + tx; i < length; i += BLOCK_SIZE) { |
| 1200 | filter_and_add_to_histogram(score[i], i); |
| 1201 | } |
| 1202 | __syncthreads(); |
| 1203 | |
| 1204 | // Stage 2: refine with 8bit radix passes |
| 1205 | #pragma unroll |
| 1206 | for (int round = 0; round < NUM_ROUNDS; ++round) { |
| 1207 | auto& s_last_remain = storage.last_remain; |
| 1208 | const auto r_idx = round % 2; |
| 1209 | |
| 1210 | const auto _raw_num_input = s_num_input[r_idx]; |
| 1211 | if constexpr (CheckOverflow) { |
| 1212 | if (_raw_num_input > SMEM_INPUT_SIZE) return false; |
| 1213 | } |
| 1214 | const auto num_input = |
| 1215 | (_raw_num_input < SMEM_INPUT_SIZE) ? _raw_num_input : SMEM_INPUT_SIZE; |
| 1216 | |
| 1217 | run_cumsum(); |
| 1218 | if (tx < RADIX && s_histogram[tx] > topk && s_histogram[tx + 1] <= topk) { |
| 1219 | s_threshold_bin_id = tx; |
| 1220 | s_num_input[r_idx ^ 1] = 0; |
| 1221 | s_last_remain = topk - s_histogram[tx + 1]; |
| 1222 | } |
| 1223 | __syncthreads(); |
| 1224 | |
| 1225 | const auto threshold = s_threshold_bin_id; |
| 1226 | topk -= s_histogram[threshold + 1]; |
| 1227 | |
| 1228 | const int offset = FIRST_SHIFT - round * 8; |
| 1229 | const bool is_last_round = (round == NUM_ROUNDS - 1); |
| 1230 | |
| 1231 | if (topk == 0) { |
| 1232 | for (int i = tx; i < num_input; i += BLOCK_SIZE) { |
| 1233 | const auto idx = s_input_idx[r_idx][i]; |
| 1234 | const auto bin = (Traits::ToOrdered(score[idx]) >> offset) & 0xFF; |
| 1235 | if (static_cast<int>(bin) > threshold) { |
| 1236 | const auto pos = atomicAdd(&s_counter, 1); |
| 1237 | s_indices[pos] = idx; |
| 1238 | } |
| 1239 | } |
| 1240 | __syncthreads(); |
| 1241 | break; |
| 1242 | } else { |
| 1243 | __syncthreads(); |
| 1244 | if (tx < RADIX + 1) s_histogram[tx] = 0; |
| 1245 | __syncthreads(); |
| 1246 | for (int i = tx; i < num_input; i += BLOCK_SIZE) { |
| 1247 | const auto idx = s_input_idx[r_idx][i]; |
| 1248 | const auto raw_input = score[idx]; |
| 1249 | const auto bin = (Traits::ToOrdered(raw_input) >> offset) & 0xFF; |
| 1250 | if (static_cast<int>(bin) > threshold) { |
| 1251 | const auto pos = atomicAdd(&s_counter, 1); |
| 1252 | s_indices[pos] = idx; |
| 1253 | } else if (static_cast<int>(bin) == threshold) { |
| 1254 | if (is_last_round) { |
| 1255 | const auto pos = atomicAdd(&s_last_remain, -1); |
| 1256 | if (pos > 0) { |
| 1257 | s_indices[top_k - pos] = idx; |
| 1258 | } |
| 1259 | } else { |
| 1260 | const auto pos = atomicAdd(&s_num_input[r_idx ^ 1], 1); |
| 1261 | if (__builtin_expect(pos < SMEM_INPUT_SIZE, 1)) { |
| 1262 | s_input_idx[r_idx ^ 1][pos] = idx; |
| 1263 | const auto bin32 = Traits::ToOrdered(raw_input); |
| 1264 | const auto sub_bin = (bin32 >> (offset - 8)) & 0xFF; |
| 1265 | atomicAdd(&s_histogram[sub_bin], 1); |
| 1266 | } |
| 1267 | } |
| 1268 | } |
| 1269 | } |
| 1270 | __syncthreads(); |
| 1271 | } |
| 1272 | } |
| 1273 | } |
| 1274 | |
| 1275 | // Output phase - mode-specific |
| 1276 | #pragma unroll 2 |
| 1277 | for (int base = tx; base < static_cast<int>(top_k); base += BLOCK_SIZE) { |
| 1278 | const int idx = s_indices[base]; |
| 1279 | dst[base] = static_cast<IdType>(idx); |
| 1280 | } |
| 1281 | return true; |
| 1282 | } |
| 1283 | |
| 1284 | /*! |
| 1285 | * \brief Filtered Top-K kernel for ragged sequences. |
| 1286 | * |
| 1287 | * \tparam DType Data type (float, half, nv_bfloat16) |
| 1288 | * \tparam IdType Index type (int32_t) |
| 1289 | * \tparam VEC_SIZE Vector size for input loads (1, 2, 4, or 8) |
| 1290 | */ |
| 1291 | template <typename DType, typename IdType, int VEC_SIZE, uint32_t MAX_K = 2048, |
| 1292 | bool UsePredicatedShortLoads = false> |
| 1293 | __global__ void __launch_bounds__(FILTERED_TOPK_BLOCK_THREADS) |
| 1294 | FilteredTopKUnifiedKernel(const DType* __restrict__ input, |
| 1295 | IdType* __restrict__ output, |
| 1296 | const IdType* __restrict__ lengths, |
| 1297 | uint32_t num_rows, uint32_t top_k, |
| 1298 | uint32_t max_len) { |
| 1299 | const uint32_t bid = blockIdx.x; |
| 1300 | if (bid >= num_rows) return; |
| 1301 | const int length = lengths ? lengths[bid] : static_cast<int>(max_len); |
| 1302 | __shared__ FilteredTopKStorage<MAX_K> storage; |
| 1303 | filtered_topk_row<DType, IdType, VEC_SIZE, MAX_K, UsePredicatedShortLoads>( |
| 1304 | input + bid * max_len, output + bid * top_k, length, top_k, storage); |
| 1305 | } |
| 1306 | |
| 1307 | // Helper to compute GCD for VEC_SIZE selection |
| 1308 | constexpr uint32_t gcd(uint32_t a, uint32_t b) { |
| 1309 | while (b != 0) { |
| 1310 | uint32_t t = b; |
| 1311 | b = a % b; |
| 1312 | a = t; |
| 1313 | } |
| 1314 | return a; |
| 1315 | } |
| 1316 | |
| 1317 | // Compute optimal VEC_SIZE based on max_len and dtype |
| 1318 | // Returns 1, 2, 4, or 8 |
| 1319 | template <typename DType> |
| 1320 | constexpr int ComputeFilteredTopKVecSize(uint32_t max_len) { |
| 1321 | constexpr int MAX_VEC = 16 / sizeof(DType); // 4 for float32, 8 for fp16/bf16 |
| 1322 | // Use GCD to find largest power-of-2 divisor |
| 1323 | const uint32_t g = gcd(max_len, static_cast<uint32_t>(MAX_VEC)); |
| 1324 | return static_cast<int>(g); |
| 1325 | } |
| 1326 | |
| 1327 | template <typename DType, typename IdType, uint32_t MAX_K = 2048> |
| 1328 | cudaError_t FilteredTopKRaggedTransform(const DType* input, |
| 1329 | IdType* output_indices, |
| 1330 | const IdType* lengths, |
| 1331 | uint32_t num_rows, uint32_t top_k_val, |
| 1332 | uint32_t max_len, |
| 1333 | cudaStream_t stream = 0) { |
| 1334 | constexpr size_t smem_size = FILTERED_TOPK_SMEM_DYNAMIC; |
| 1335 | constexpr int MAX_VEC = 16 / sizeof(DType); |
| 1336 | |
| 1337 | dim3 grid(num_rows); |
| 1338 | dim3 block(FILTERED_TOPK_BLOCK_THREADS); |
| 1339 | void* args[] = {&input, &output_indices, &lengths, |
| 1340 | &num_rows, &top_k_val, &max_len}; |
| 1341 | |
| 1342 | const int vec_size = ComputeFilteredTopKVecSize<DType>(max_len); |
| 1343 | |
| 1344 | #define DISPATCH_VEC_SIZE(VS) \ |
| 1345 | if (vec_size == VS) { \ |
| 1346 | auto kernel = \ |
| 1347 | FilteredTopKUnifiedKernel<DType, IdType, VS, MAX_K, (VS != MAX_VEC)>; \ |
| 1348 | FLASHINFER_CUDA_CALL(cudaFuncSetAttribute( \ |
| 1349 | kernel, cudaFuncAttributeMaxDynamicSharedMemorySize, smem_size)); \ |
| 1350 | FLASHINFER_CUDA_CALL(cudaLaunchKernel((void*)kernel, grid, block, args, \ |
| 1351 | smem_size, stream)); \ |
| 1352 | return cudaSuccess; \ |
| 1353 | } |
| 1354 | |
| 1355 | DISPATCH_VEC_SIZE(1) |
| 1356 | DISPATCH_VEC_SIZE(2) |
| 1357 | DISPATCH_VEC_SIZE(4) |
| 1358 | if constexpr (MAX_VEC >= 8) { |
| 1359 | DISPATCH_VEC_SIZE(8) |
| 1360 | } |
| 1361 | #undef DISPATCH_VEC_SIZE |
| 1362 | |
| 1363 | return cudaSuccess; |
| 1364 | } |
| 1365 | |
| 1366 | } // namespace filtered_topk |
| 1367 | |
| 1368 | template <typename DType, typename IdType, uint32_t MAX_K = 2048> |
| 1369 | cudaError_t FilteredTopKRaggedTransform(const DType* input, |
| 1370 | IdType* output_indices, |
| 1371 | const IdType* lengths, |
| 1372 | uint32_t num_rows, uint32_t top_k_val, |
| 1373 | uint32_t max_len, |
| 1374 | cudaStream_t stream = 0) { |
| 1375 | return filtered_topk::FilteredTopKRaggedTransform<DType, IdType, MAX_K>( |
| 1376 | input, output_indices, lengths, num_rows, top_k_val, max_len, stream); |
| 1377 | } |
| 1378 | |
| 1379 | } // namespace vllm |
| 1380 | |
| 1381 | #endif // PERSISTENT_TOPK_CUH_ |