| 1 | // Adapted from |
| 2 | // https://github.com/sgl-project/sglang/tree/main/sgl-kernel/csrc/cpu: |
| 3 | // `biased_topk_kernel_impl`/`biased_topk_cpu` (flat, non-grouped DeepSeek-V4 |
| 4 | // biased top-k) and `hash_topk_kernel_impl`/`hash_topk_cpu` (hash-routed |
| 5 | // layers, expert IDs from a precomputed `tid2eid` lookup table), plus their |
| 6 | // `sigmoid`/`softmax`/`apply_bias`/`sqrtsoftplus` scoring helpers. See |
| 7 | // `biased_topk_cpu`/`hash_topk_cpu` below for the one deliberate deviation |
| 8 | // from the upstream kernel (routed_scaling_factor handling when no shared |
| 9 | // expert is fused into routing, vLLM's CPU DeepSeek-V4 calling convention). |
| 10 | // |
| 11 | // clang-format off |
| 12 | |
| 13 | #include <algorithm> |
| 14 | #include <array> |
| 15 | #include <cmath> |
| 16 | #include <limits> |
| 17 | #include <type_traits> |
| 18 | #include <utility> |
| 19 | |
| 20 | #include "common.h" |
| 21 | #include "vec.h" |
| 22 | |
| 23 | namespace { |
| 24 | |
| 25 | template <typename scalar_t, int SIZE, std::enable_if_t<!std::is_same_v<scalar_t, float>, int> = 0> |
| 26 | inline void sigmoid(float* __restrict__ out, const scalar_t* __restrict__ input) { |
| 27 | using bVec = at::vec::Vectorized<scalar_t>; |
| 28 | using fVec = at::vec::Vectorized<float>; |
| 29 | |
| 30 | const fVec one = fVec(1.f); |
| 31 | |
| 32 | constexpr int kVecSize = bVec::size(); |
| 33 | for (int d = 0; d < SIZE; d += kVecSize) { |
| 34 | bVec x_bvec = bVec::loadu(input + d); |
| 35 | fVec x_fvec0, x_fvec1; |
| 36 | std::tie(x_fvec0, x_fvec1) = at::vec::convert_to_float(x_bvec); |
| 37 | |
| 38 | x_fvec0 = one / (one + x_fvec0.neg().exp_u20()); |
| 39 | x_fvec1 = one / (one + x_fvec1.neg().exp_u20()); |
| 40 | |
| 41 | x_fvec0.store(out + d); |
| 42 | x_fvec1.store(out + d + fVec::size()); |
| 43 | } |
| 44 | } |
| 45 | |
| 46 | template <typename scalar_t, int SIZE, std::enable_if_t<std::is_same_v<scalar_t, float>, int> = 0> |
| 47 | inline void sigmoid(float* __restrict__ out, const float* __restrict__ input) { |
| 48 | using fVec = at::vec::Vectorized<float>; |
| 49 | const fVec one = fVec(1.f); |
| 50 | constexpr int kVecSize = fVec::size(); |
| 51 | for (int d = 0; d < SIZE; d += kVecSize) { |
| 52 | fVec in_fvec = fVec::loadu(input + d); |
| 53 | in_fvec = one / (one + in_fvec.neg().exp_u20()); |
| 54 | in_fvec.store(out + d); |
| 55 | } |
| 56 | } |
| 57 | |
| 58 | template <typename scalar_t, int SIZE, std::enable_if_t<!std::is_same_v<scalar_t, float>, int> = 0> |
| 59 | inline void softmax(float* __restrict__ out, const scalar_t* __restrict__ input) { |
| 60 | using bVec = at::vec::Vectorized<scalar_t>; |
| 61 | using fVec = at::vec::Vectorized<float>; |
| 62 | |
| 63 | constexpr int kVecSize = bVec::size(); |
| 64 | |
| 65 | // step 1: get max |
| 66 | fVec max_fvec = fVec(-std::numeric_limits<float>::infinity()); |
| 67 | for (int d = 0; d < SIZE; d += kVecSize) { |
| 68 | bVec x_bvec = bVec::loadu(input + d); |
| 69 | fVec x_fvec0, x_fvec1; |
| 70 | std::tie(x_fvec0, x_fvec1) = at::vec::convert_to_float(x_bvec); |
| 71 | |
| 72 | max_fvec = at::vec::maximum(max_fvec, x_fvec0); |
| 73 | max_fvec = at::vec::maximum(max_fvec, x_fvec1); |
| 74 | x_fvec0.store(out + d); |
| 75 | x_fvec1.store(out + d + fVec::size()); |
| 76 | } |
| 77 | float max_val = vec_reduce_max(max_fvec); |
| 78 | max_fvec = fVec(max_val); |
| 79 | |
| 80 | // step 2: sum of (x - max).exp() |
| 81 | fVec sum_fvec = fVec(float(0)); |
| 82 | for (int d = 0; d < SIZE; d += fVec::size()) { |
| 83 | fVec x_fvec = (fVec::loadu(out + d) - max_fvec).exp_u20(); |
| 84 | sum_fvec += x_fvec; |
| 85 | x_fvec.store(out + d); |
| 86 | } |
| 87 | float sum_val = vec_reduce_sum(sum_fvec); |
| 88 | |
| 89 | // step 3: x * (1 / sum) |
| 90 | sum_fvec = fVec(1.f / sum_val); |
| 91 | for (int d = 0; d < SIZE; d += fVec::size()) { |
| 92 | fVec out_fvec = fVec::loadu(out + d) * sum_fvec; |
| 93 | out_fvec.store(out + d); |
| 94 | } |
| 95 | } |
| 96 | |
| 97 | template <typename scalar_t, int SIZE, std::enable_if_t<std::is_same_v<scalar_t, float>, int> = 0> |
| 98 | inline void softmax(float* __restrict__ out, const float* __restrict__ input) { |
| 99 | using fVec = at::vec::Vectorized<float>; |
| 100 | |
| 101 | constexpr int kVecSize = fVec::size(); |
| 102 | |
| 103 | // step 1: get max |
| 104 | fVec max_fvec = fVec(-std::numeric_limits<float>::infinity()); |
| 105 | for (int d = 0; d < SIZE; d += kVecSize) { |
| 106 | fVec x_fvec = fVec::loadu(input + d); |
| 107 | max_fvec = at::vec::maximum(max_fvec, x_fvec); |
| 108 | x_fvec.store(out + d); |
| 109 | } |
| 110 | float max_val = vec_reduce_max(max_fvec); |
| 111 | max_fvec = fVec(max_val); |
| 112 | |
| 113 | // step 2: sum of (x - max).exp() |
| 114 | fVec sum_fvec = fVec(float(0)); |
| 115 | for (int d = 0; d < SIZE; d += kVecSize) { |
| 116 | fVec x_fvec = (fVec::loadu(out + d) - max_fvec).exp_u20(); |
| 117 | sum_fvec += x_fvec; |
| 118 | x_fvec.store(out + d); |
| 119 | } |
| 120 | float sum_val = vec_reduce_sum(sum_fvec); |
| 121 | |
| 122 | // step 3: x * (1 / sum) |
| 123 | sum_fvec = fVec(1.f / sum_val); |
| 124 | for (int d = 0; d < SIZE; d += kVecSize) { |
| 125 | fVec out_fvec = fVec::loadu(out + d) * sum_fvec; |
| 126 | out_fvec.store(out + d); |
| 127 | } |
| 128 | } |
| 129 | |
| 130 | template <typename param_t, int SIZE> |
| 131 | inline void |
| 132 | apply_bias(float* __restrict__ scores2, const float* __restrict__ scores, const param_t* __restrict__ bias) { |
| 133 | using fVec = at::vec::Vectorized<float>; |
| 134 | auto vec_size = fVec::size() * 2; |
| 135 | int d = 0; |
| 136 | for (; d <= SIZE - vec_size; d += vec_size) { |
| 137 | fVec bias0, bias1, x0, x1; |
| 138 | std::tie(bias0, bias1) = load_float_vec2(bias + d); |
| 139 | std::tie(x0, x1) = load_float_vec2(scores + d); |
| 140 | x0 = x0 + bias0; |
| 141 | x1 = x1 + bias1; |
| 142 | x0.store(scores2 + d); |
| 143 | x1.store(scores2 + d + fVec::size()); |
| 144 | } |
| 145 | for (; d < SIZE; d++) { |
| 146 | scores2[d] = scores[d] + (float)bias[d]; |
| 147 | } |
| 148 | } |
| 149 | |
| 150 | // sqrtsoftplus: sqrt(softplus(x)) = sqrt(log(1 + exp(x))) |
| 151 | // For numerical stability: when x > threshold, softplus(x) ≈ x |
| 152 | // When x < -threshold, softplus(x) ≈ exp(x), so sqrt(softplus(x)) ≈ sqrt(exp(x)) = exp(x/2) |
| 153 | template <typename scalar_t, int SIZE, std::enable_if_t<!std::is_same_v<scalar_t, float>, int> = 0> |
| 154 | inline void sqrtsoftplus(float* __restrict__ out, const scalar_t* __restrict__ input) { |
| 155 | using bVec = at::vec::Vectorized<scalar_t>; |
| 156 | using fVec = at::vec::Vectorized<float>; |
| 157 | |
| 158 | const fVec one = fVec(1.f); |
| 159 | const fVec half = fVec(0.5f); |
| 160 | const fVec threshold = fVec(20.f); |
| 161 | const fVec neg_threshold = fVec(-15.f); // below this, 1+exp(x) loses precision in float32 |
| 162 | |
| 163 | constexpr int kVecSize = bVec::size(); |
| 164 | for (int d = 0; d < SIZE; d += kVecSize) { |
| 165 | bVec x_bvec = bVec::loadu(input + d); |
| 166 | fVec x_fvec0, x_fvec1; |
| 167 | std::tie(x_fvec0, x_fvec1) = at::vec::convert_to_float(x_bvec); |
| 168 | |
| 169 | // default: softplus(x) = log(1 + exp(x)) |
| 170 | fVec sp0 = (one + x_fvec0.exp_u20()).log(); |
| 171 | fVec sp1 = (one + x_fvec1.exp_u20()).log(); |
| 172 | // x > 20: softplus(x) ≈ x |
| 173 | sp0 = fVec::blendv(sp0, x_fvec0, x_fvec0 > threshold); |
| 174 | sp1 = fVec::blendv(sp1, x_fvec1, x_fvec1 > threshold); |
| 175 | // x < -15: softplus(x) ≈ exp(x), sqrt(exp(x)) = exp(x/2) |
| 176 | fVec exp_half0 = (x_fvec0 * half).exp_u20(); |
| 177 | fVec exp_half1 = (x_fvec1 * half).exp_u20(); |
| 178 | sp0 = fVec::blendv(sp0.sqrt(), exp_half0, x_fvec0 < neg_threshold); |
| 179 | sp1 = fVec::blendv(sp1.sqrt(), exp_half1, x_fvec1 < neg_threshold); |
| 180 | |
| 181 | sp0.store(out + d); |
| 182 | sp1.store(out + d + fVec::size()); |
| 183 | } |
| 184 | } |
| 185 | |
| 186 | template <typename scalar_t, int SIZE, std::enable_if_t<std::is_same_v<scalar_t, float>, int> = 0> |
| 187 | inline void sqrtsoftplus(float* __restrict__ out, const float* __restrict__ input) { |
| 188 | using fVec = at::vec::Vectorized<float>; |
| 189 | const fVec one = fVec(1.f); |
| 190 | const fVec half = fVec(0.5f); |
| 191 | const fVec threshold = fVec(20.f); |
| 192 | const fVec neg_threshold = fVec(-15.f); |
| 193 | constexpr int kVecSize = fVec::size(); |
| 194 | for (int d = 0; d < SIZE; d += kVecSize) { |
| 195 | fVec x = fVec::loadu(input + d); |
| 196 | fVec sp = (one + x.exp_u20()).log(); |
| 197 | sp = fVec::blendv(sp, x, x > threshold); |
| 198 | fVec exp_half = (x * half).exp_u20(); |
| 199 | sp = fVec::blendv(sp.sqrt(), exp_half, x < neg_threshold); |
| 200 | sp.store(out + d); |
| 201 | } |
| 202 | } |
| 203 | |
| 204 | // biased_topk: flat (non-grouped) biased topk for DeepSeek V4 |
| 205 | // scoring_func: 0 = sigmoid, 1 = sqrtsoftplus |
| 206 | template <typename scalar_t, typename param_t, int NUM_EXPERTS, int TOPK> |
| 207 | void biased_topk_kernel_impl( |
| 208 | float* __restrict__ topk_weights, |
| 209 | int32_t* __restrict__ topk_ids, |
| 210 | const scalar_t* __restrict__ gating_output, |
| 211 | const param_t* __restrict__ bias, |
| 212 | int64_t num_tokens, |
| 213 | int64_t topk, |
| 214 | bool renormalize, |
| 215 | int scoring_func, |
| 216 | int64_t num_fused_shared_experts, |
| 217 | float routed_scaling_factor, |
| 218 | bool apply_routed_scaling_factor_on_output) { |
| 219 | at::parallel_for(0, num_tokens, 0, [&](int64_t begin, int64_t end) { |
| 220 | alignas(64) float scores[NUM_EXPERTS]; |
| 221 | alignas(64) float scores_biased[NUM_EXPERTS]; |
| 222 | |
| 223 | using elem_t = std::pair<float, int32_t>; |
| 224 | std::array<elem_t, NUM_EXPERTS> queue; |
| 225 | |
| 226 | // simple RNG for fused shared expert random ID |
| 227 | uint64_t rng_state = begin * 6364136223846793005ULL + 1442695040888963407ULL; |
| 228 | |
| 229 | for (int64_t i = begin; i < end; ++i) { |
| 230 | // compute scores |
| 231 | if (scoring_func == 0) { |
| 232 | sigmoid<scalar_t, NUM_EXPERTS>(scores, gating_output + i * NUM_EXPERTS); |
| 233 | } else { |
| 234 | sqrtsoftplus<scalar_t, NUM_EXPERTS>(scores, gating_output + i * NUM_EXPERTS); |
| 235 | } |
| 236 | |
| 237 | // add bias for selection |
| 238 | apply_bias<param_t, NUM_EXPERTS>(scores_biased, scores, bias); |
| 239 | |
| 240 | // build queue and partial sort to find top-k |
| 241 | for (int64_t e = 0; e < NUM_EXPERTS; ++e) { |
| 242 | queue[e] = {scores_biased[e], static_cast<int32_t>(e)}; |
| 243 | } |
| 244 | |
| 245 | std::partial_sort( |
| 246 | queue.begin(), queue.begin() + topk, queue.end(), [](const elem_t& x, const elem_t& y) -> bool { |
| 247 | return x.first > y.first; |
| 248 | }); |
| 249 | |
| 250 | // gather original scores (without bias) as weights |
| 251 | for (int64_t j = 0; j < topk; ++j) { |
| 252 | int32_t idx = queue[j].second; |
| 253 | topk_ids[i * topk + j] = idx; |
| 254 | topk_weights[i * topk + j] = scores[idx]; |
| 255 | } |
| 256 | |
| 257 | // handle fused shared experts |
| 258 | if (num_fused_shared_experts > 0) { |
| 259 | // replace last slot with random shared expert ID |
| 260 | rng_state = rng_state * 6364136223846793005ULL + 1442695040888963407ULL; |
| 261 | int32_t shared_id = |
| 262 | NUM_EXPERTS + static_cast<int32_t>((rng_state >> 33) % static_cast<uint64_t>(num_fused_shared_experts)); |
| 263 | topk_ids[i * topk + topk - 1] = shared_id; |
| 264 | |
| 265 | // shared expert weight = sum of routed weights / scaling_factor |
| 266 | if (routed_scaling_factor != 0.0f) { |
| 267 | float routed_sum = 0.f; |
| 268 | for (int64_t j = 0; j < topk - 1; ++j) { |
| 269 | routed_sum += topk_weights[i * topk + j]; |
| 270 | } |
| 271 | topk_weights[i * topk + topk - 1] = routed_sum / routed_scaling_factor; |
| 272 | } |
| 273 | } |
| 274 | |
| 275 | // renormalize |
| 276 | if (renormalize) { |
| 277 | float sum = 0.f; |
| 278 | int64_t norm_end = (num_fused_shared_experts == 0) ? topk : topk - 1; |
| 279 | for (int64_t j = 0; j < norm_end; ++j) { |
| 280 | sum += topk_weights[i * topk + j]; |
| 281 | } |
| 282 | float scale = 1.f / sum; |
| 283 | if (apply_routed_scaling_factor_on_output) { |
| 284 | scale *= routed_scaling_factor; |
| 285 | } |
| 286 | for (int64_t j = 0; j < norm_end; ++j) { |
| 287 | topk_weights[i * topk + j] *= scale; |
| 288 | } |
| 289 | // also scale the shared expert weight if present |
| 290 | if (num_fused_shared_experts > 0) { |
| 291 | topk_weights[i * topk + topk - 1] *= scale; |
| 292 | } |
| 293 | } |
| 294 | } |
| 295 | }); |
| 296 | } |
| 297 | |
| 298 | // hash_topk: expert IDs come from a precomputed lookup table tid2eid[input_ids] |
| 299 | // scoring_func: 0 = softmax, 1 = sigmoid, 2 = sqrtsoftplus |
| 300 | template <typename scalar_t, int NUM_EXPERTS, int TOPK> |
| 301 | void hash_topk_kernel_impl( |
| 302 | float* __restrict__ topk_weights, |
| 303 | int32_t* __restrict__ topk_ids, |
| 304 | const scalar_t* __restrict__ gating_output, |
| 305 | const int32_t* __restrict__ tid2eid, // [num_tokens, routed_topk] |
| 306 | int64_t num_tokens, |
| 307 | int scoring_func, |
| 308 | int64_t num_fused_shared_experts, |
| 309 | int64_t num_experts, |
| 310 | float routed_scaling_factor, |
| 311 | int64_t topk) { |
| 312 | const int64_t routed_topk = topk - num_fused_shared_experts; |
| 313 | const bool need_renormalize = (scoring_func != 0); // renormalize for non-softmax |
| 314 | |
| 315 | at::parallel_for(0, num_tokens, 0, [&](int64_t begin, int64_t end) { |
| 316 | alignas(64) float scores[NUM_EXPERTS]; |
| 317 | |
| 318 | // simple RNG for fused shared expert random ID |
| 319 | uint64_t rng_state = begin * 6364136223846793005ULL + 1442695040888963407ULL; |
| 320 | |
| 321 | for (int64_t i = begin; i < end; ++i) { |
| 322 | // compute scores over all experts |
| 323 | if (scoring_func == 0) { |
| 324 | softmax<scalar_t, NUM_EXPERTS>(scores, gating_output + i * NUM_EXPERTS); |
| 325 | } else if (scoring_func == 1) { |
| 326 | sigmoid<scalar_t, NUM_EXPERTS>(scores, gating_output + i * NUM_EXPERTS); |
| 327 | } else { |
| 328 | sqrtsoftplus<scalar_t, NUM_EXPERTS>(scores, gating_output + i * NUM_EXPERTS); |
| 329 | } |
| 330 | |
| 331 | // gather expert IDs from lookup table |
| 332 | const int32_t* eid_row = tid2eid + i * routed_topk; |
| 333 | for (int64_t j = 0; j < routed_topk; ++j) { |
| 334 | int32_t eid = eid_row[j]; |
| 335 | topk_ids[i * topk + j] = eid; |
| 336 | topk_weights[i * topk + j] = scores[eid]; |
| 337 | } |
| 338 | |
| 339 | // renormalize routed weights (for non-softmax scoring). Unlike the |
| 340 | // upstream kernel, always fold routed_scaling_factor in here rather |
| 341 | // than only applying it to the fused-shared-expert slot below: vLLM's |
| 342 | // DeepSeek-V4 CPU hash-routed layers never fuse a shared expert into |
| 343 | // this routing table (num_fused_shared_experts is always 0 on this |
| 344 | // call path), so this is the only place the scaling factor can be |
| 345 | // applied. Matches the eager reference (`sqrtsoftplus_bias_topk` in |
| 346 | // cpu_moe.py) and the CUDA kernel (`dsv4HashTopkSoftplusSqrt`). |
| 347 | if (need_renormalize) { |
| 348 | float sum = 0.f; |
| 349 | for (int64_t j = 0; j < routed_topk; ++j) { |
| 350 | sum += topk_weights[i * topk + j]; |
| 351 | } |
| 352 | if (sum > 0.f) { |
| 353 | float scale = 1.f / sum; |
| 354 | if (num_fused_shared_experts == 0) { |
| 355 | scale *= routed_scaling_factor; |
| 356 | } |
| 357 | for (int64_t j = 0; j < routed_topk; ++j) { |
| 358 | topk_weights[i * topk + j] *= scale; |
| 359 | } |
| 360 | } |
| 361 | } |
| 362 | |
| 363 | // handle fused shared expert |
| 364 | if (num_fused_shared_experts > 0) { |
| 365 | // random shared expert ID in [num_experts, num_experts + num_fused_shared_experts) |
| 366 | rng_state = rng_state * 6364136223846793005ULL + 1442695040888963407ULL; |
| 367 | int32_t shared_id = |
| 368 | num_experts + static_cast<int32_t>((rng_state >> 33) % static_cast<uint64_t>(num_fused_shared_experts)); |
| 369 | topk_ids[i * topk + topk - 1] = shared_id; |
| 370 | |
| 371 | // shared expert weight = sum of routed weights / scaling_factor |
| 372 | float routed_sum = 0.f; |
| 373 | for (int64_t j = 0; j < routed_topk; ++j) { |
| 374 | routed_sum += topk_weights[i * topk + j]; |
| 375 | } |
| 376 | topk_weights[i * topk + topk - 1] = routed_sum / routed_scaling_factor; |
| 377 | } |
| 378 | } |
| 379 | }); |
| 380 | } |
| 381 | |
| 382 | #define LAUNCH_HASH_TOPK_KERNEL(NE, NTOPK) \ |
| 383 | hash_topk_kernel_impl<scalar_t, NE, NTOPK>( \ |
| 384 | topk_weights.data_ptr<float>(), \ |
| 385 | topk_ids.data_ptr<int32_t>(), \ |
| 386 | gating_output.data_ptr<scalar_t>(), \ |
| 387 | tid2eid_flat.data_ptr<int32_t>(), \ |
| 388 | num_tokens, \ |
| 389 | scoring_func_id, \ |
| 390 | num_fused_shared_experts, \ |
| 391 | num_experts, \ |
| 392 | routed_scaling_factor_value, \ |
| 393 | topk); |
| 394 | |
| 395 | #define LAUNCH_BIASED_TOPK_KERNEL(NE, NTOPK) \ |
| 396 | biased_topk_kernel_impl<scalar_t, param_t, NE, NTOPK>( \ |
| 397 | topk_weights.data_ptr<float>(), \ |
| 398 | topk_ids.data_ptr<int32_t>(), \ |
| 399 | gating_output.data_ptr<scalar_t>(), \ |
| 400 | correction_bias.data_ptr<param_t>(), \ |
| 401 | num_tokens, \ |
| 402 | topk, \ |
| 403 | renormalize, \ |
| 404 | scoring_func_id, \ |
| 405 | num_fused_shared_experts, \ |
| 406 | routed_scaling_factor_value, \ |
| 407 | apply_routed_scaling_factor_on_output); |
| 408 | |
| 409 | } // namespace |
| 410 | |
| 411 | // biased topk for DeepSeek V4 (flat, non-grouped) |
| 412 | std::tuple<at::Tensor, at::Tensor> biased_topk_cpu( |
| 413 | at::Tensor& hidden_states, |
| 414 | at::Tensor& gating_output, |
| 415 | at::Tensor& correction_bias, |
| 416 | int64_t topk, |
| 417 | bool renormalize, |
| 418 | std::string scoring_func, |
| 419 | int64_t num_fused_shared_experts, |
| 420 | std::optional<double> routed_scaling_factor, |
| 421 | bool apply_routed_scaling_factor_on_output) { |
| 422 | CHECK_INPUT(gating_output); |
| 423 | CHECK_INPUT(correction_bias); |
| 424 | |
| 425 | const auto st = gating_output.scalar_type(); |
| 426 | int64_t num_tokens = hidden_states.size(0); |
| 427 | int64_t num_experts = gating_output.size(1); |
| 428 | TORCH_CHECK(gating_output.size(0) == num_tokens, "Number of tokens mismatch"); |
| 429 | TORCH_CHECK(correction_bias.numel() == num_experts, "Bias shape mismatch"); |
| 430 | |
| 431 | int scoring_func_id = 0; |
| 432 | if (scoring_func == "sigmoid") { |
| 433 | scoring_func_id = 0; |
| 434 | } else if (scoring_func == "sqrtsoftplus") { |
| 435 | scoring_func_id = 1; |
| 436 | } else { |
| 437 | TORCH_CHECK(false, "Unsupported scoring_func: ", scoring_func); |
| 438 | } |
| 439 | |
| 440 | float routed_scaling_factor_value = routed_scaling_factor.has_value() ? routed_scaling_factor.value() : 0.0f; |
| 441 | |
| 442 | at::Tensor topk_weights = at::empty({num_tokens, topk}, hidden_states.options().dtype(at::kFloat)); |
| 443 | at::Tensor topk_ids = at::empty({num_tokens, topk}, hidden_states.options().dtype(at::kInt)); |
| 444 | |
| 445 | // The actual routed topk (excluding fused shared experts) |
| 446 | int64_t routed_topk = topk - num_fused_shared_experts; |
| 447 | |
| 448 | CPU_DISPATCH_FLOATING_TYPES_EXT(st, correction_bias.scalar_type(), "biased_topk_kernel", [&] { |
| 449 | // dispatch on num_experts and routed_topk. For DeepSeek V4: 256 experts, |
| 450 | // routed_topk=6 (num_fused_shared_experts is always 0 on this CPU call |
| 451 | // path -- see cpu_moe.py). Anything outside this table isn't a real |
| 452 | // model shape, so fail loudly instead of silently degrading to a slow |
| 453 | // generic path. |
| 454 | switch (num_experts) { |
| 455 | case 64: |
| 456 | switch (routed_topk) { |
| 457 | case 6: |
| 458 | LAUNCH_BIASED_TOPK_KERNEL(64, 6); |
| 459 | break; |
| 460 | case 8: |
| 461 | LAUNCH_BIASED_TOPK_KERNEL(64, 8); |
| 462 | break; |
| 463 | default: |
| 464 | TORCH_CHECK(false, "biased_topk_cpu: unsupported routed_topk ", routed_topk, |
| 465 | " for num_experts=64"); |
| 466 | } |
| 467 | break; |
| 468 | case 128: |
| 469 | switch (routed_topk) { |
| 470 | case 6: |
| 471 | LAUNCH_BIASED_TOPK_KERNEL(128, 6); |
| 472 | break; |
| 473 | case 8: |
| 474 | LAUNCH_BIASED_TOPK_KERNEL(128, 8); |
| 475 | break; |
| 476 | default: |
| 477 | TORCH_CHECK(false, "biased_topk_cpu: unsupported routed_topk ", routed_topk, |
| 478 | " for num_experts=128"); |
| 479 | } |
| 480 | break; |
| 481 | case 256: |
| 482 | switch (routed_topk) { |
| 483 | case 6: |
| 484 | LAUNCH_BIASED_TOPK_KERNEL(256, 6); |
| 485 | break; |
| 486 | case 8: |
| 487 | LAUNCH_BIASED_TOPK_KERNEL(256, 8); |
| 488 | break; |
| 489 | case 9: |
| 490 | LAUNCH_BIASED_TOPK_KERNEL(256, 9); |
| 491 | break; |
| 492 | default: |
| 493 | TORCH_CHECK(false, "biased_topk_cpu: unsupported routed_topk ", routed_topk, |
| 494 | " for num_experts=256"); |
| 495 | } |
| 496 | break; |
| 497 | case 384: |
| 498 | switch (routed_topk) { |
| 499 | case 6: |
| 500 | LAUNCH_BIASED_TOPK_KERNEL(384, 6); |
| 501 | break; |
| 502 | case 8: |
| 503 | LAUNCH_BIASED_TOPK_KERNEL(384, 8); |
| 504 | break; |
| 505 | default: |
| 506 | TORCH_CHECK(false, "biased_topk_cpu: unsupported routed_topk ", routed_topk, |
| 507 | " for num_experts=384"); |
| 508 | } |
| 509 | break; |
| 510 | default: |
| 511 | TORCH_CHECK(false, "biased_topk_cpu: unsupported num_experts ", num_experts); |
| 512 | } |
| 513 | }); |
| 514 | return std::make_tuple(topk_weights, topk_ids); |
| 515 | } |
| 516 | |
| 517 | // hash topk for DeepSeek V4 (expert IDs from precomputed lookup table) |
| 518 | std::tuple<at::Tensor, at::Tensor> hash_topk_cpu( |
| 519 | at::Tensor& gating_output, |
| 520 | at::Tensor& tid2eid, |
| 521 | int64_t topk, |
| 522 | std::string scoring_func, |
| 523 | int64_t num_fused_shared_experts, |
| 524 | int64_t num_experts, |
| 525 | double routed_scaling_factor) { |
| 526 | CHECK_INPUT(gating_output); |
| 527 | CHECK_INPUT(tid2eid); |
| 528 | |
| 529 | const auto st = gating_output.scalar_type(); |
| 530 | int64_t num_tokens = gating_output.size(0); |
| 531 | int64_t num_experts_gating = gating_output.size(1); |
| 532 | int64_t routed_topk = topk - num_fused_shared_experts; |
| 533 | |
| 534 | TORCH_CHECK(tid2eid.size(0) == num_tokens, "tid2eid row count must match num_tokens"); |
| 535 | TORCH_CHECK(tid2eid.size(1) == routed_topk, "tid2eid column count must match routed_topk"); |
| 536 | TORCH_CHECK(tid2eid.scalar_type() == at::kInt, "tid2eid must be int32"); |
| 537 | TORCH_CHECK(num_experts_gating == num_experts, "num_experts mismatch"); |
| 538 | |
| 539 | int scoring_func_id = 0; |
| 540 | if (scoring_func == "softmax") { |
| 541 | scoring_func_id = 0; |
| 542 | } else if (scoring_func == "sigmoid") { |
| 543 | scoring_func_id = 1; |
| 544 | } else if (scoring_func == "sqrtsoftplus") { |
| 545 | scoring_func_id = 2; |
| 546 | } else { |
| 547 | TORCH_CHECK(false, "Unsupported scoring_func: ", scoring_func); |
| 548 | } |
| 549 | |
| 550 | float routed_scaling_factor_value = static_cast<float>(routed_scaling_factor); |
| 551 | |
| 552 | at::Tensor topk_weights = at::empty({num_tokens, topk}, gating_output.options().dtype(at::kFloat)); |
| 553 | at::Tensor topk_ids = at::empty({num_tokens, topk}, gating_output.options().dtype(at::kInt)); |
| 554 | |
| 555 | // tid2eid is [num_tokens, routed_topk], already indexed by input_ids in |
| 556 | // Python; CHECK_INPUT above already guarantees it is contiguous. |
| 557 | const at::Tensor& tid2eid_flat = tid2eid; |
| 558 | |
| 559 | // Dispatch for bf16, fp16, and float32 |
| 560 | // Note: cannot use AT_DISPATCH_FLOATING_TYPES_AND2 since it includes double which lacks convert_to_float |
| 561 | auto dispatch_fn = [&]<typename scalar_t>() { |
| 562 | // For DeepSeek V4: 256 experts, routed_topk in {6,7,8} depending on |
| 563 | // layer (num_fused_shared_experts is always 0 on this CPU call path -- |
| 564 | // see cpu_moe.py). Anything outside this table isn't a real model |
| 565 | // shape, so fail loudly instead of silently degrading to a slow |
| 566 | // generic path. |
| 567 | switch (num_experts) { |
| 568 | case 64: |
| 569 | switch (routed_topk) { |
| 570 | case 6: |
| 571 | LAUNCH_HASH_TOPK_KERNEL(64, 6); |
| 572 | break; |
| 573 | case 7: |
| 574 | LAUNCH_HASH_TOPK_KERNEL(64, 7); |
| 575 | break; |
| 576 | case 8: |
| 577 | LAUNCH_HASH_TOPK_KERNEL(64, 8); |
| 578 | break; |
| 579 | default: |
| 580 | TORCH_CHECK(false, "hash_topk_cpu: unsupported routed_topk ", routed_topk, |
| 581 | " for num_experts=64"); |
| 582 | } |
| 583 | break; |
| 584 | case 128: |
| 585 | switch (routed_topk) { |
| 586 | case 6: |
| 587 | LAUNCH_HASH_TOPK_KERNEL(128, 6); |
| 588 | break; |
| 589 | case 7: |
| 590 | LAUNCH_HASH_TOPK_KERNEL(128, 7); |
| 591 | break; |
| 592 | case 8: |
| 593 | LAUNCH_HASH_TOPK_KERNEL(128, 8); |
| 594 | break; |
| 595 | default: |
| 596 | TORCH_CHECK(false, "hash_topk_cpu: unsupported routed_topk ", routed_topk, |
| 597 | " for num_experts=128"); |
| 598 | } |
| 599 | break; |
| 600 | case 256: |
| 601 | switch (routed_topk) { |
| 602 | case 6: |
| 603 | LAUNCH_HASH_TOPK_KERNEL(256, 6); |
| 604 | break; |
| 605 | case 7: |
| 606 | LAUNCH_HASH_TOPK_KERNEL(256, 7); |
| 607 | break; |
| 608 | case 8: |
| 609 | LAUNCH_HASH_TOPK_KERNEL(256, 8); |
| 610 | break; |
| 611 | default: |
| 612 | TORCH_CHECK(false, "hash_topk_cpu: unsupported routed_topk ", routed_topk, |
| 613 | " for num_experts=256"); |
| 614 | } |
| 615 | break; |
| 616 | case 384: |
| 617 | switch (routed_topk) { |
| 618 | case 6: |
| 619 | LAUNCH_HASH_TOPK_KERNEL(384, 6); |
| 620 | break; |
| 621 | case 7: |
| 622 | LAUNCH_HASH_TOPK_KERNEL(384, 7); |
| 623 | break; |
| 624 | case 8: |
| 625 | LAUNCH_HASH_TOPK_KERNEL(384, 8); |
| 626 | break; |
| 627 | default: |
| 628 | TORCH_CHECK(false, "hash_topk_cpu: unsupported routed_topk ", routed_topk, |
| 629 | " for num_experts=384"); |
| 630 | } |
| 631 | break; |
| 632 | default: |
| 633 | TORCH_CHECK(false, "hash_topk_cpu: unsupported num_experts ", num_experts); |
| 634 | } |
| 635 | }; |
| 636 | |
| 637 | if (st == at::ScalarType::BFloat16) { |
| 638 | dispatch_fn.template operator()<at::BFloat16>(); |
| 639 | } else if (st == at::ScalarType::Half) { |
| 640 | dispatch_fn.template operator()<at::Half>(); |
| 641 | } else if (st == at::ScalarType::Float) { |
| 642 | dispatch_fn.template operator()<float>(); |
| 643 | } else { |
| 644 | TORCH_CHECK(false, "Unsupported dtype for hash_topk_cpu: ", st); |
| 645 | } |
| 646 | |
| 647 | return std::make_tuple(topk_weights, topk_ids); |
| 648 | } |
| 649 | |
| 650 | // Adapted from |
| 651 | // https://github.com/sgl-project/sglang/tree/main/sgl-kernel/csrc/cpu: |
| 652 | // `topk_transform_512_cpu`/`topk_transform_512_cpu_kernel_impl`. Near- |
| 653 | // verbatim port for the sparse indexer's DECODE path: hardcodes top-512 |
| 654 | // (the only `index_topk` DeepSeek-V4-Flash uses) and streams a min-heap over |
| 655 | // the tail of the row instead of sorting the whole row. |
| 656 | // |
| 657 | // Differences from upstream: |
| 658 | // - `out_page_indices` is still computed (required by the kernel contract) |
| 659 | // but vLLM's caller only consumes `out_raw_indices` (local/compressed- |
| 660 | // context positions -- `DeepseekV4CPUAttention.forward_mqa` resolves those |
| 661 | // to physical slots itself via `map_local_to_global_slots_cpu`). |
| 662 | // - `page_size` was already a runtime argument upstream, unlike |
| 663 | // `fp8_paged_mqa_logits_cpu` (paged_mqa_logits.cpp), which needed one. |
| 664 | namespace { |
| 665 | |
| 666 | constexpr int64_t kC4Topk = 512; |
| 667 | |
| 668 | template <typename T> |
| 669 | inline int64_t load_index_value(const T* __restrict__ ptr, int64_t idx) { |
| 670 | return static_cast<int64_t>(ptr[idx]); |
| 671 | } |
| 672 | |
| 673 | struct TopKTransformElem { |
| 674 | float score; |
| 675 | int32_t index; |
| 676 | }; |
| 677 | |
| 678 | struct TopKTransformMinHeapCmp { |
| 679 | bool operator()(const TopKTransformElem& lhs, const TopKTransformElem& rhs) const { |
| 680 | if (lhs.score == rhs.score) { |
| 681 | return lhs.index > rhs.index; |
| 682 | } |
| 683 | return lhs.score > rhs.score; |
| 684 | } |
| 685 | }; |
| 686 | |
| 687 | template <typename seq_t, typename page_t> |
| 688 | void topk_transform_512_cpu_kernel_impl( |
| 689 | const float* __restrict__ scores, |
| 690 | const seq_t* __restrict__ seq_lens, |
| 691 | const page_t* __restrict__ page_tables, |
| 692 | int32_t* __restrict__ out_page_indices, |
| 693 | int32_t* __restrict__ out_raw_indices, |
| 694 | int64_t batch_size, |
| 695 | int64_t max_seq_len, |
| 696 | int64_t page_table_stride, |
| 697 | int64_t out_stride, |
| 698 | int64_t page_size) { |
| 699 | TORCH_CHECK(page_size > 0, "page_size must be positive"); |
| 700 | TORCH_CHECK((page_size & (page_size - 1)) == 0, "page_size must be a power of 2"); |
| 701 | const int page_bits = page_size > 1 ? static_cast<int>(std::log2(static_cast<double>(page_size))) : 0; |
| 702 | const int64_t page_mask = page_size - 1; |
| 703 | |
| 704 | at::parallel_for(0, batch_size, 0, [&](int64_t begin, int64_t end) { |
| 705 | std::array<TopKTransformElem, kC4Topk> heap; |
| 706 | |
| 707 | for (int64_t b = begin; b < end; ++b) { |
| 708 | const float* __restrict__ scores_row = scores + b * max_seq_len; |
| 709 | const page_t* __restrict__ page_table_row = page_tables + b * page_table_stride; |
| 710 | int32_t* __restrict__ out_page_row = out_page_indices + b * out_stride; |
| 711 | int32_t* __restrict__ out_raw_row = out_raw_indices == nullptr ? nullptr : out_raw_indices + b * out_stride; |
| 712 | |
| 713 | int64_t seq_len = load_index_value(seq_lens, b); |
| 714 | seq_len = std::max<int64_t>(0, std::min<int64_t>(seq_len, max_seq_len)); |
| 715 | const int64_t valid_topk = std::min<int64_t>(seq_len, kC4Topk); |
| 716 | |
| 717 | auto store_slot = [&](int64_t slot, int32_t raw_index) { |
| 718 | if (raw_index < 0) { |
| 719 | out_page_row[slot] = -1; |
| 720 | if (out_raw_row != nullptr) { |
| 721 | out_raw_row[slot] = -1; |
| 722 | } |
| 723 | return; |
| 724 | } |
| 725 | |
| 726 | const int64_t page_idx = static_cast<int64_t>(raw_index) >> page_bits; |
| 727 | const int64_t offset_in_page = static_cast<int64_t>(raw_index) & page_mask; |
| 728 | const int64_t physical_page = load_index_value(page_table_row, page_idx); |
| 729 | out_page_row[slot] = static_cast<int32_t>((physical_page << page_bits) | offset_in_page); |
| 730 | if (out_raw_row != nullptr) { |
| 731 | out_raw_row[slot] = raw_index; |
| 732 | } |
| 733 | }; |
| 734 | |
| 735 | if (seq_len <= kC4Topk) { |
| 736 | for (int64_t i = 0; i < valid_topk; ++i) { |
| 737 | store_slot(i, static_cast<int32_t>(i)); |
| 738 | } |
| 739 | for (int64_t i = valid_topk; i < kC4Topk; ++i) { |
| 740 | store_slot(i, -1); |
| 741 | } |
| 742 | continue; |
| 743 | } |
| 744 | |
| 745 | for (int64_t i = 0; i < kC4Topk; ++i) { |
| 746 | heap[i] = {scores_row[i], static_cast<int32_t>(i)}; |
| 747 | } |
| 748 | std::make_heap(heap.begin(), heap.end(), TopKTransformMinHeapCmp()); |
| 749 | |
| 750 | for (int64_t i = kC4Topk; i < seq_len; ++i) { |
| 751 | const float score = scores_row[i]; |
| 752 | const TopKTransformElem& current_min = heap.front(); |
| 753 | if (score > current_min.score || |
| 754 | (score == current_min.score && static_cast<int32_t>(i) < current_min.index)) { |
| 755 | std::pop_heap(heap.begin(), heap.end(), TopKTransformMinHeapCmp()); |
| 756 | heap.back() = {score, static_cast<int32_t>(i)}; |
| 757 | std::push_heap(heap.begin(), heap.end(), TopKTransformMinHeapCmp()); |
| 758 | } |
| 759 | } |
| 760 | |
| 761 | for (int64_t i = 0; i < kC4Topk; ++i) { |
| 762 | store_slot(i, heap[i].index); |
| 763 | } |
| 764 | } |
| 765 | }); |
| 766 | } |
| 767 | |
| 768 | } // namespace |
| 769 | |
| 770 | void topk_transform_512_cpu( |
| 771 | at::Tensor& scores, |
| 772 | at::Tensor& seq_lens, |
| 773 | at::Tensor& page_tables, |
| 774 | at::Tensor& out_page_indices, |
| 775 | int64_t page_size, |
| 776 | const std::optional<at::Tensor>& out_raw_indices) { |
| 777 | CHECK_INPUT(scores); |
| 778 | CHECK_INPUT(seq_lens); |
| 779 | CHECK_INPUT(page_tables); |
| 780 | CHECK_INPUT(out_page_indices); |
| 781 | |
| 782 | TORCH_CHECK(scores.scalar_type() == at::kFloat, "scores must be float32"); |
| 783 | TORCH_CHECK(out_page_indices.scalar_type() == at::kInt, "out_page_indices must be int32"); |
| 784 | TORCH_CHECK(scores.dim() == 2, "scores must be a 2D tensor"); |
| 785 | TORCH_CHECK(seq_lens.dim() == 1, "seq_lens must be a 1D tensor"); |
| 786 | TORCH_CHECK(page_tables.dim() == 2, "page_tables must be a 2D tensor"); |
| 787 | TORCH_CHECK(out_page_indices.dim() == 2, "out_page_indices must be a 2D tensor"); |
| 788 | |
| 789 | const int64_t batch_size = scores.size(0); |
| 790 | const int64_t max_seq_len = scores.size(1); |
| 791 | TORCH_CHECK(seq_lens.size(0) == batch_size, "seq_lens row count must match scores"); |
| 792 | TORCH_CHECK(page_tables.size(0) == batch_size, "page_tables row count must match scores"); |
| 793 | TORCH_CHECK(out_page_indices.size(0) == batch_size, "out_page_indices row count must match scores"); |
| 794 | TORCH_CHECK(out_page_indices.size(1) >= kC4Topk, "out_page_indices must have at least 512 columns"); |
| 795 | |
| 796 | int32_t* raw_ptr = nullptr; |
| 797 | if (out_raw_indices.has_value()) { |
| 798 | at::Tensor raw = out_raw_indices.value(); |
| 799 | CHECK_INPUT(raw); |
| 800 | TORCH_CHECK(raw.scalar_type() == at::kInt, "out_raw_indices must be int32"); |
| 801 | TORCH_CHECK(raw.dim() == 2, "out_raw_indices must be a 2D tensor"); |
| 802 | TORCH_CHECK(raw.sizes() == out_page_indices.sizes(), "out_raw_indices shape must match out_page_indices"); |
| 803 | raw_ptr = raw.data_ptr<int32_t>(); |
| 804 | } |
| 805 | |
| 806 | if (batch_size == 0) { |
| 807 | return; |
| 808 | } |
| 809 | |
| 810 | auto launch_with_seq_type = [&](auto seq_zero) { |
| 811 | using seq_t = decltype(seq_zero); |
| 812 | if (page_tables.scalar_type() == at::kInt) { |
| 813 | topk_transform_512_cpu_kernel_impl<seq_t, int32_t>( |
| 814 | scores.data_ptr<float>(), seq_lens.data_ptr<seq_t>(), page_tables.data_ptr<int32_t>(), |
| 815 | out_page_indices.data_ptr<int32_t>(), raw_ptr, batch_size, max_seq_len, page_tables.stride(0), |
| 816 | out_page_indices.stride(0), page_size); |
| 817 | } else if (page_tables.scalar_type() == at::kLong) { |
| 818 | topk_transform_512_cpu_kernel_impl<seq_t, int64_t>( |
| 819 | scores.data_ptr<float>(), seq_lens.data_ptr<seq_t>(), page_tables.data_ptr<int64_t>(), |
| 820 | out_page_indices.data_ptr<int32_t>(), raw_ptr, batch_size, max_seq_len, page_tables.stride(0), |
| 821 | out_page_indices.stride(0), page_size); |
| 822 | } else { |
| 823 | TORCH_CHECK(false, "page_tables must be int32 or int64"); |
| 824 | } |
| 825 | }; |
| 826 | |
| 827 | if (seq_lens.scalar_type() == at::kInt) { |
| 828 | launch_with_seq_type(int32_t{0}); |
| 829 | } else if (seq_lens.scalar_type() == at::kLong) { |
| 830 | launch_with_seq_type(int64_t{0}); |
| 831 | } else { |
| 832 | TORCH_CHECK(false, "seq_lens must be int32 or int64"); |
| 833 | } |
| 834 | } |