main
cu 901 lines 39.5 KB
Raw
1 /*
2 * Adapted from https://github.com/NVIDIA/TensorRT-LLM/blob/v0.7.1/cpp/tensorrt_llm/kernels/mixtureOfExperts/moe_kernels.cu
3 * Copyright (c) 2024, The vLLM team.
4 * SPDX-FileCopyrightText: Copyright (c) 1993-2023 NVIDIA CORPORATION & AFFILIATES. All rights reserved.
5 * SPDX-License-Identifier: Apache-2.0
6 *
7 * Licensed under the Apache License, Version 2.0 (the "License");
8 * you may not use this file except in compliance with the License.
9 * You may obtain a copy of the License at
10 *
11 * http://www.apache.org/licenses/LICENSE-2.0
12 *
13 * Unless required by applicable law or agreed to in writing, software
14 * distributed under the License is distributed on an "AS IS" BASIS,
15 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16 * See the License for the specific language governing permissions and
17 * limitations under the License.
18 */
19 #include <type_traits>
20
21 #include <cuda_runtime.h>
22 #include <torch/csrc/stable/accelerator.h>
23 #include <torch/csrc/stable/tensor.h>
24 #include <torch/headeronly/core/ScalarType.h>
25 #include <torch/headeronly/util/Exception.h>
26
27 #include "../../cuda_compat.h"
28 #include "../cub_helpers.h"
29 #include "libtorch_stable/torch_utils.h"
30
31 #ifndef USE_ROCM
32 #include <cuda_bf16.h>
33 #include <cuda_fp16.h>
34 #else
35 #include <hip/hip_bf16.h>
36 #include <hip/hip_fp16.h>
37 typedef __hip_bfloat16 __nv_bfloat16;
38 typedef __hip_bfloat162 __nv_bfloat162;
39 #endif
40
41 #define MAX(a, b) ((a) > (b) ? (a) : (b))
42 #define MIN(a, b) ((a) < (b) ? (a) : (b))
43
44 namespace vllm {
45 namespace moe {
46
47 /// Aligned array type
48 template <
49 typename T,
50 /// Number of elements in the array
51 int N,
52 /// Alignment requirement in bytes
53 int Alignment = sizeof(T) * N
54 >
55 struct alignas(Alignment) AlignedArray {
56 T data[N];
57 };
58
59 template <typename T>
60 __device__ __forceinline__ float toFloat(T value) {
61 if constexpr (std::is_same_v<T, float>) {
62 return value;
63 } else if constexpr (std::is_same_v<T, __nv_bfloat16>) {
64 return __bfloat162float(value);
65 } else if constexpr (std::is_same_v<T, __half>) {
66 return __half2float(value);
67 }
68 }
69
70 // Scoring function enums
71 enum ScoringFunc {
72 SCORING_SOFTMAX = 0, // apply softmax
73 SCORING_SIGMOID = 1 // apply sigmoid
74 };
75
76 // ====================== Softmax things ===============================
77 // We have our own implementation of softmax here so we can support transposing the output
78 // in the softmax kernel when we extend this module to support expert-choice routing.
79 template <int TPB, typename InputType>
80 __launch_bounds__(TPB) __global__
81 void moeSoftmax(const InputType* input, const bool* finished, float* output, const int num_cols)
82 {
83 using BlockReduce = cub::BlockReduce<float, TPB>;
84 __shared__ typename BlockReduce::TempStorage tmpStorage;
85
86 __shared__ float normalizing_factor;
87 __shared__ float float_max;
88
89 const int thread_row_offset = blockIdx.x * num_cols;
90
91 float threadData(-FLT_MAX);
92
93 // Don't touch finished rows.
94 if ((finished != nullptr) && finished[blockIdx.x])
95 {
96 return;
97 }
98
99 for (int ii = threadIdx.x; ii < num_cols; ii += TPB)
100 {
101 const int idx = thread_row_offset + ii;
102 const float val = toFloat(input[idx]);
103 threadData = max(val, threadData);
104 }
105
106 const float maxElem = BlockReduce(tmpStorage).Reduce(threadData, CubMaxOp());
107 if (threadIdx.x == 0)
108 {
109 float_max = maxElem;
110 }
111 __syncthreads();
112
113 threadData = 0;
114
115 for (int ii = threadIdx.x; ii < num_cols; ii += TPB)
116 {
117 const int idx = thread_row_offset + ii;
118 const float val = toFloat(input[idx]);
119 threadData += expf(val - float_max);
120 }
121
122 const auto Z = BlockReduce(tmpStorage).Reduce(threadData, CubAddOp());
123
124 if (threadIdx.x == 0)
125 {
126 normalizing_factor = 1.f / Z;
127 }
128 __syncthreads();
129
130 for (int ii = threadIdx.x; ii < num_cols; ii += TPB)
131 {
132 const int idx = thread_row_offset + ii;
133 const float val = toFloat(input[idx]);
134 float softmax_val = expf(val - float_max) * normalizing_factor;
135 // Clamp NaN/Inf to 0 to prevent duplicate expert IDs downstream.
136 if (isnan(softmax_val) || isinf(softmax_val)) softmax_val = 0.f;
137 output[idx] = softmax_val;
138 }
139 }
140
141 template <int TPB, typename InputType>
142 __launch_bounds__(TPB) __global__
143 void moeSigmoid(const InputType* input, const bool* finished, float* output, const int num_cols)
144 {
145 const int thread_row_offset = blockIdx.x * num_cols;
146
147 // Don't touch finished rows.
148 if ((finished != nullptr) && finished[blockIdx.x])
149 {
150 return;
151 }
152
153 for (int ii = threadIdx.x; ii < num_cols; ii += TPB)
154 {
155 const int idx = thread_row_offset + ii;
156 const float val = toFloat(input[idx]);
157 float sigmoid_val = 1.0f / (1.0f + __expf(-val));
158 // Clamp NaN/Inf to 0 to prevent duplicate expert IDs downstream.
159 if (isnan(sigmoid_val) || isinf(sigmoid_val)) sigmoid_val = 0.f;
160 output[idx] = sigmoid_val;
161 }
162 }
163
164 template <int TPB, typename IndType>
165 __launch_bounds__(TPB) __global__ void moeTopK(
166 const float* inputs_after_softmax,
167 const bool* finished,
168 float* output,
169 IndType* indices,
170 int* source_rows,
171 const int num_experts,
172 const int k,
173 const int start_expert,
174 const int end_expert,
175 const bool renormalize,
176 const float* bias,
177 const double routed_scaling_factor,
178 const bool* is_padding)
179 {
180
181 using cub_kvp = cub::KeyValuePair<int, float>;
182 using BlockReduce = cub::BlockReduce<cub_kvp, TPB>;
183 __shared__ typename BlockReduce::TempStorage tmpStorage;
184
185 cub_kvp thread_kvp;
186 cub::ArgMax arg_max;
187
188 const int num_rows = gridDim.x;
189 const int block_row = blockIdx.x;
190
191 const bool row_is_active = finished ? !finished[block_row] : true;
192 const int thread_read_offset = blockIdx.x * num_experts;
193 float selected_sum = 0.f;
194 for (int k_idx = 0; k_idx < k; ++k_idx)
195 {
196 thread_kvp.key = 0;
197 thread_kvp.value = -1.f; // This is OK because inputs are probabilities
198
199 cub_kvp inp_kvp;
200 for (int expert = threadIdx.x; expert < num_experts; expert += TPB)
201 {
202 const int idx = thread_read_offset + expert;
203 inp_kvp.key = expert;
204
205 // Apply correction bias if provided
206 if (bias != nullptr) {
207 inp_kvp.value = inputs_after_softmax[idx] + bias[expert];
208 } else {
209 inp_kvp.value = inputs_after_softmax[idx];
210 }
211
212 for (int prior_k = 0; prior_k < k_idx; ++prior_k)
213 {
214 const int prior_winning_expert = indices[k * block_row + prior_k];
215
216 if (prior_winning_expert == expert)
217 {
218 inp_kvp = thread_kvp;
219 }
220 }
221
222 thread_kvp = arg_max(inp_kvp, thread_kvp);
223 }
224
225 const cub_kvp result_kvp = BlockReduce(tmpStorage).Reduce(thread_kvp, arg_max);
226 if (threadIdx.x == 0)
227 {
228 // Ignore experts the node isn't responsible for with expert parallelism
229 const int expert = result_kvp.key;
230 const bool node_uses_expert = expert >= start_expert && expert < end_expert;
231 const bool should_process_row = row_is_active && node_uses_expert;
232 const bool is_pad_row = is_padding != nullptr && is_padding[block_row];
233
234 const int idx = k * block_row + k_idx;
235 // Return the unbiased scores for output weights
236 output[idx] = inputs_after_softmax[thread_read_offset + expert];
237 indices[idx] = is_pad_row ? static_cast<IndType>(-1)
238 : (should_process_row ? (expert - start_expert) : num_experts);
239 assert(is_pad_row || indices[idx] >= 0);
240 source_rows[idx] = k_idx * num_rows + block_row;
241 if (renormalize) {
242 selected_sum += inputs_after_softmax[thread_read_offset + expert];
243 }
244 }
245 __syncthreads();
246 }
247
248 // Apply renormalization and routed scaling factor to final weights.
249 if (threadIdx.x == 0) {
250 float scale = static_cast<float>(routed_scaling_factor);
251 if (renormalize) {
252 const float denom = selected_sum > 0.f ? selected_sum : 1.f;
253 scale /= denom;
254 }
255 for (int k_idx = 0; k_idx < k; ++k_idx) {
256 const int idx = k * block_row + k_idx;
257 output[idx] = output[idx] * scale;
258 }
259 }
260 }
261
262 // ====================== TopK softmax things ===============================
263
264 /*
265 A Top-K gating softmax written to exploit when the number of experts in the MoE layers
266 are a small power of 2. This allows us to cleanly share the rows among the threads in
267 a single warp and eliminate communication between warps (so no need to use shared mem).
268
269 It fuses the softmax, max and argmax into a single kernel.
270
271 Limitations:
272 1) This implementation is optimized for when the number of experts is a small power of 2.
273 Additionally it also supports when number of experts is multiple of 64 which is still
274 faster than the computing softmax and topK separately (only tested on CUDA yet).
275 2) This implementation assumes k is small, but will work for any k.
276 */
277
278 template <int VPT, int NUM_EXPERTS, int WARPS_PER_CTA, int BYTES_PER_LDG, int WARP_SIZE_PARAM, typename IndType,
279 typename InputType = float, ScoringFunc SF>
280 __launch_bounds__(WARPS_PER_CTA* WARP_SIZE_PARAM) __global__
281 void topkGating(const InputType* input, const bool* finished, float* output, const int num_rows, IndType* indices,
282 int* source_rows, const int k, const int start_expert, const int end_expert, const bool renormalize,
283 const float* bias, const double routed_scaling_factor, const bool* is_padding)
284 {
285 static_assert(std::is_same_v<InputType, float> || std::is_same_v<InputType, __nv_bfloat16> ||
286 std::is_same_v<InputType, __half>,
287 "InputType must be float, __nv_bfloat16, or __half");
288
289 // We begin by enforcing compile time assertions and setting up compile time constants.
290 static_assert(BYTES_PER_LDG == (BYTES_PER_LDG & -BYTES_PER_LDG), "BYTES_PER_LDG must be power of 2");
291 static_assert(BYTES_PER_LDG <= 16, "BYTES_PER_LDG must be leq 16");
292
293 // Number of bytes each thread pulls in per load
294 static constexpr int ELTS_PER_LDG = BYTES_PER_LDG / sizeof(InputType);
295 static constexpr int ELTS_PER_ROW = NUM_EXPERTS;
296 static constexpr int THREADS_PER_ROW = ELTS_PER_ROW / VPT;
297 static constexpr int LDG_PER_THREAD = VPT / ELTS_PER_LDG;
298
299 if constexpr (std::is_same_v<InputType, __nv_bfloat16> || std::is_same_v<InputType, __half>) {
300 static_assert(ELTS_PER_LDG == 1 || ELTS_PER_LDG % 2 == 0,
301 "ELTS_PER_LDG must be 1 or even for 16-bit conversion");
302 }
303
304 // Restrictions based on previous section.
305 static_assert(VPT % ELTS_PER_LDG == 0, "The elements per thread must be a multiple of the elements per ldg");
306 static_assert(WARP_SIZE_PARAM % THREADS_PER_ROW == 0, "The threads per row must cleanly divide the threads per warp");
307 static_assert(THREADS_PER_ROW == (THREADS_PER_ROW & -THREADS_PER_ROW), "THREADS_PER_ROW must be power of 2");
308 static_assert(THREADS_PER_ROW <= WARP_SIZE_PARAM, "THREADS_PER_ROW can be at most warp size");
309
310 // We have NUM_EXPERTS elements per row. We specialize for small #experts
311 static constexpr int ELTS_PER_WARP = WARP_SIZE_PARAM * VPT;
312 static constexpr int ROWS_PER_WARP = ELTS_PER_WARP / ELTS_PER_ROW;
313 static constexpr int ROWS_PER_CTA = WARPS_PER_CTA * ROWS_PER_WARP;
314
315 // Restrictions for previous section.
316 static_assert(ELTS_PER_WARP % ELTS_PER_ROW == 0, "The elts per row must cleanly divide the total elt per warp");
317
318 // ===================== From this point, we finally start computing run-time variables. ========================
319
320 // Compute CTA and warp rows. We pack multiple rows into a single warp, and a block contains WARPS_PER_CTA warps.
321 // This, each block processes a chunk of rows. We start by computing the start row for each block.
322 const int cta_base_row = blockIdx.x * ROWS_PER_CTA;
323
324 // Now, using the base row per thread block, we compute the base row per warp.
325 const int warp_base_row = cta_base_row + threadIdx.y * ROWS_PER_WARP;
326
327 // The threads in a warp are split into sub-groups that will work on a row.
328 // We compute row offset for each thread sub-group
329 const int thread_row_in_warp = threadIdx.x / THREADS_PER_ROW;
330 const int thread_row = warp_base_row + thread_row_in_warp;
331
332 // Threads with indices out of bounds should early exit here.
333 if (thread_row >= num_rows)
334 {
335 return;
336 }
337 const bool row_is_active = finished ? !finished[thread_row] : true;
338
339 // We finally start setting up the read pointers for each thread. First, each thread jumps to the start of the
340 // row it will read.
341 const InputType* thread_row_ptr = input + thread_row * ELTS_PER_ROW;
342
343 // Now, we compute the group each thread belong to in order to determine the first column to start loads.
344 const int thread_group_idx = threadIdx.x % THREADS_PER_ROW;
345 const int first_elt_read_by_thread = thread_group_idx * ELTS_PER_LDG;
346 const InputType* thread_read_ptr = thread_row_ptr + first_elt_read_by_thread;
347
348 // Finally, we pull in the data from global mem
349 float row_chunk[VPT];
350
351 // NOTE(zhuhaoran): dispatch different input types loading, BF16/FP16 convert to float
352 if constexpr (std::is_same_v<InputType, float>) {
353 using VecType = AlignedArray<float, ELTS_PER_LDG>;
354 VecType* row_chunk_vec_ptr = reinterpret_cast<VecType*>(&row_chunk);
355 const VecType* vec_thread_read_ptr = reinterpret_cast<const VecType*>(thread_read_ptr);
356 #pragma unroll
357 for (int ii = 0; ii < LDG_PER_THREAD; ++ii) {
358 row_chunk_vec_ptr[ii] = vec_thread_read_ptr[ii * THREADS_PER_ROW];
359 }
360 } else if constexpr (std::is_same_v<InputType, __nv_bfloat16>) {
361 if constexpr (ELTS_PER_LDG >= 2) {
362 using VecType = AlignedArray<__nv_bfloat16, ELTS_PER_LDG>;
363 float2* row_chunk_f2 = reinterpret_cast<float2*>(row_chunk);
364 const VecType* vec_thread_read_ptr = reinterpret_cast<const VecType*>(thread_read_ptr);
365 #pragma unroll
366 for (int ii = 0; ii < LDG_PER_THREAD; ++ii) {
367 VecType vec = vec_thread_read_ptr[ii * THREADS_PER_ROW];
368 int base_idx_f2 = ii * ELTS_PER_LDG / 2;
369 #pragma unroll
370 for (int jj = 0; jj < ELTS_PER_LDG / 2; ++jj) {
371 row_chunk_f2[base_idx_f2 + jj] = __bfloat1622float2(
372 *reinterpret_cast<const __nv_bfloat162*>(vec.data + jj * 2)
373 );
374 }
375 }
376 } else { // ELTS_PER_LDG == 1
377 #pragma unroll
378 for (int ii = 0; ii < LDG_PER_THREAD; ++ii) {
379 const __nv_bfloat16* scalar_ptr = thread_read_ptr + ii * THREADS_PER_ROW;
380 row_chunk[ii] = __bfloat162float(*scalar_ptr);
381 }
382 }
383 } else if constexpr (std::is_same_v<InputType, __half>) {
384 if constexpr (ELTS_PER_LDG >= 2) {
385 using VecType = AlignedArray<__half, ELTS_PER_LDG>;
386 float2* row_chunk_f2 = reinterpret_cast<float2*>(row_chunk);
387 const VecType* vec_thread_read_ptr = reinterpret_cast<const VecType*>(thread_read_ptr);
388 #pragma unroll
389 for (int ii = 0; ii < LDG_PER_THREAD; ++ii) {
390 VecType vec = vec_thread_read_ptr[ii * THREADS_PER_ROW];
391 int base_idx_f2 = ii * ELTS_PER_LDG / 2;
392 #pragma unroll
393 for (int jj = 0; jj < ELTS_PER_LDG / 2; ++jj) {
394 row_chunk_f2[base_idx_f2 + jj] = __half22float2(
395 *reinterpret_cast<const __half2*>(vec.data + jj * 2)
396 );
397 }
398 }
399 } else { // ELTS_PER_LDG == 1
400 #pragma unroll
401 for (int ii = 0; ii < LDG_PER_THREAD; ++ii) {
402 const __half* scalar_ptr = thread_read_ptr + ii * THREADS_PER_ROW;
403 row_chunk[ii] = __half2float(*scalar_ptr);
404 }
405 }
406 }
407
408 if constexpr (SF == SCORING_SOFTMAX) {
409 // First, we perform a max reduce within the thread.
410 float thread_max = row_chunk[0];
411 #pragma unroll
412 for (int ii = 1; ii < VPT; ++ii) {
413 thread_max = max(thread_max, row_chunk[ii]);
414 }
415
416 // Now, we find the max within the thread group and distribute among the threads. We use a butterfly reduce.
417 #pragma unroll
418 for (int mask = THREADS_PER_ROW / 2; mask > 0; mask /= 2)
419 {
420 thread_max = max(thread_max, VLLM_SHFL_XOR_SYNC_WIDTH(thread_max, mask, THREADS_PER_ROW));
421 }
422
423 // From this point, thread max in all the threads have the max within the row.
424 // Now, we subtract the max from each element in the thread and take the exp. We also compute the thread local sum.
425 float row_sum = 0;
426 #pragma unroll
427 for (int ii = 0; ii < VPT; ++ii)
428 {
429 row_chunk[ii] = expf(row_chunk[ii] - thread_max);
430 row_sum += row_chunk[ii];
431 }
432
433 // Now, we perform the sum reduce within each thread group. Similar to the max reduce, we use a bufferfly pattern.
434 #pragma unroll
435 for (int mask = THREADS_PER_ROW / 2; mask > 0; mask /= 2)
436 {
437 row_sum += VLLM_SHFL_XOR_SYNC_WIDTH(row_sum, mask, THREADS_PER_ROW);
438 }
439
440 // From this point, all threads have the max and the sum for their rows in the thread_max and thread_sum variables
441 // respectively. Finally, we can scale the rows for the softmax. Technically, for top-k gating we don't need to
442 // compute the entire softmax row. We can likely look at the maxes and only compute for the top-k values in the row.
443 // However, this kernel will likely not be a bottle neck and it seems better to closer match torch and find the
444 // argmax after computing the softmax.
445 const float reciprocal_row_sum = 1.f / row_sum;
446
447 #pragma unroll
448 for (int ii = 0; ii < VPT; ++ii)
449 {
450 row_chunk[ii] = row_chunk[ii] * reciprocal_row_sum;
451 }
452 } else if constexpr (SF == SCORING_SIGMOID) {
453 #pragma unroll
454 for (int ii = 0; ii < VPT; ++ii)
455 {
456 row_chunk[ii] = 1.0f / (1.0f + __expf(-row_chunk[ii]));
457 }
458 }
459
460 // Fix: clamp NaN/Inf values to 0 to prevent duplicate expert IDs.
461 // NaN gating (from degenerate hidden states in CUDA graph padding) causes
462 // softmax to produce all-NaN, which makes the argmax loop always pick
463 // expert 0 for every top-k slot, producing duplicate expert IDs that
464 // crash FlashInfer's three-step MoE sort.
465 // With 0s, the argmax uses index tie-breaking to pick [0,1,2,...,k-1].
466 #pragma unroll
467 for (int ii = 0; ii < VPT; ++ii) {
468 if (isnan(row_chunk[ii]) || isinf(row_chunk[ii])) {
469 row_chunk[ii] = 0.f;
470 }
471 }
472
473 static constexpr int COLS_PER_GROUP_LDG = ELTS_PER_LDG * THREADS_PER_ROW;
474
475 // If bias is not null, use biased value for selection
476 float row_chunk_for_choice[VPT];
477 // Apply correction bias
478 if (bias != nullptr) {
479 #pragma unroll
480 for (int ldg = 0; ldg < LDG_PER_THREAD; ++ldg) {
481 #pragma unroll
482 for (int ii = 0; ii < ELTS_PER_LDG; ++ii) {
483 const int expert = first_elt_read_by_thread + ldg * COLS_PER_GROUP_LDG + ii;
484 float bias_val = expert < NUM_EXPERTS ? bias[expert] : 0.0f;
485 row_chunk_for_choice[ldg * ELTS_PER_LDG + ii] = row_chunk[ldg * ELTS_PER_LDG + ii] + bias_val;
486 }
487 }
488 } else {
489 #pragma unroll
490 for (int ii = 0; ii < VPT; ++ii) {
491 row_chunk_for_choice[ii] = row_chunk[ii];
492 }
493 }
494
495 // Now, row_chunk contains the softmax / sigmoid of the row chunk. Now, I want to find the topk elements in each row, along
496 // with the max index.
497 int start_col = first_elt_read_by_thread;
498
499 float selected_sum = 0.f;
500 for (int k_idx = 0; k_idx < k; ++k_idx)
501 {
502 // First, each thread does the local argmax
503 float max_val_for_choice = row_chunk_for_choice[0];
504 float max_val = row_chunk[0];
505 int expert = start_col;
506 #pragma unroll
507 for (int ldg = 0, col = start_col; ldg < LDG_PER_THREAD; ++ldg, col += COLS_PER_GROUP_LDG)
508 {
509 #pragma unroll
510 for (int ii = 0; ii < ELTS_PER_LDG; ++ii)
511 {
512 float val_for_choice = row_chunk_for_choice[ldg * ELTS_PER_LDG + ii];
513 float val = row_chunk[ldg * ELTS_PER_LDG + ii];
514
515 // No check on the experts here since columns with the smallest index are processed first and only
516 // updated if > (not >=)
517 if (val_for_choice > max_val_for_choice)
518 {
519 max_val_for_choice = val_for_choice;
520 max_val = val;
521 expert = col + ii;
522 }
523 }
524 }
525
526 // Now, we perform the argmax reduce. We use the butterfly pattern so threads reach consensus about the max.
527 // This will be useful for K > 1 so that the threads can agree on "who" had the max value. That thread can
528 // then blank out their max with -inf and the warp can run more iterations...
529 #pragma unroll
530 for (int mask = THREADS_PER_ROW / 2; mask > 0; mask /= 2)
531 {
532 float other_max_for_choice = VLLM_SHFL_XOR_SYNC_WIDTH(max_val_for_choice, mask, THREADS_PER_ROW);
533 float other_max = VLLM_SHFL_XOR_SYNC_WIDTH(max_val, mask, THREADS_PER_ROW);
534 int other_expert = VLLM_SHFL_XOR_SYNC_WIDTH(expert, mask, THREADS_PER_ROW);
535
536 // We want lower indices to "win" in every thread so we break ties this way
537 if (other_max_for_choice > max_val_for_choice || (other_max_for_choice == max_val_for_choice && other_expert < expert))
538 {
539 max_val_for_choice = other_max_for_choice;
540 max_val = other_max;
541 expert = other_expert;
542 }
543 }
544
545 // Write the max for this k iteration to global memory.
546 if (thread_group_idx == 0)
547 {
548 // Add a guard to ignore experts not included by this node
549 const bool node_uses_expert = expert >= start_expert && expert < end_expert;
550 const bool should_process_row = row_is_active && node_uses_expert;
551 const bool is_pad_row = is_padding != nullptr && is_padding[thread_row];
552
553 // The lead thread from each sub-group will write out the final results to global memory. (This will be a
554 // single) thread per row of the input/output matrices.
555 const int idx = k * thread_row + k_idx;
556 output[idx] = max_val;
557 indices[idx] = is_pad_row ? static_cast<IndType>(-1)
558 : (should_process_row ? (expert - start_expert) : NUM_EXPERTS);
559 source_rows[idx] = k_idx * num_rows + thread_row;
560 if (renormalize) {
561 selected_sum += max_val;
562 }
563 }
564
565 // Finally, we clear the value in the thread with the current max if there is another iteration to run.
566 if (k_idx + 1 < k)
567 {
568 const int ldg_group_for_expert = expert / COLS_PER_GROUP_LDG;
569 const int thread_to_clear_in_group = (expert / ELTS_PER_LDG) % THREADS_PER_ROW;
570
571 // Only the thread in the group which produced the max will reset the "winning" value to -inf.
572 if (thread_group_idx == thread_to_clear_in_group)
573 {
574 const int offset_for_expert = expert % ELTS_PER_LDG;
575 // Safe to set to any negative value since row_chunk values must be between 0 and 1.
576 row_chunk_for_choice[ldg_group_for_expert * ELTS_PER_LDG + offset_for_expert] = -10000.f;
577 }
578 }
579 }
580
581 // Apply renormalization and routed scaling factor to final weights.
582 if (thread_group_idx == 0) {
583 float scale = static_cast<float>(routed_scaling_factor);
584 if (renormalize) {
585 const float denom = selected_sum > 0.f ? selected_sum : 1.f;
586 scale /= denom;
587 }
588 for (int k_idx = 0; k_idx < k; ++k_idx) {
589 const int idx = k * thread_row + k_idx;
590 output[idx] = output[idx] * scale;
591 }
592 }
593 }
594
595 namespace detail
596 {
597 // Constructs some constants needed to partition the work across threads at compile time.
598 template <int EXPERTS, int BYTES_PER_LDG, int WARP_SIZE_PARAM, typename InputType>
599 struct TopkConstants
600 {
601 static constexpr int ELTS_PER_LDG = BYTES_PER_LDG / sizeof(InputType);
602 static_assert(EXPERTS / (ELTS_PER_LDG * WARP_SIZE_PARAM) == 0 || EXPERTS % (ELTS_PER_LDG * WARP_SIZE_PARAM) == 0, "");
603 static constexpr int VECs_PER_THREAD = MAX(1, EXPERTS / (ELTS_PER_LDG * WARP_SIZE_PARAM));
604 static constexpr int VPT = VECs_PER_THREAD * ELTS_PER_LDG;
605 static constexpr int THREADS_PER_ROW = EXPERTS / VPT;
606 static const int ROWS_PER_WARP = WARP_SIZE_PARAM / THREADS_PER_ROW;
607 };
608 } // namespace detail
609
610 template <int EXPERTS, int WARPS_PER_TB, int WARP_SIZE_PARAM, int MAX_BYTES_PER_LDG, typename IndType, typename InputType, ScoringFunc SF>
611 void topkGatingLauncherHelper(const InputType* input, const bool* finished, float* output, IndType* indices,
612 int* source_row, const int num_rows, const int k, const int start_expert, const int end_expert, const bool renormalize,
613 const float* bias, const double routed_scaling_factor, cudaStream_t stream, const bool* is_padding)
614 {
615 static constexpr int BYTES_PER_LDG = MIN(MAX_BYTES_PER_LDG, sizeof(InputType) * EXPERTS);
616 using Constants = detail::TopkConstants<EXPERTS, BYTES_PER_LDG, WARP_SIZE_PARAM, InputType>;
617 static constexpr int VPT = Constants::VPT;
618 static constexpr int ROWS_PER_WARP = Constants::ROWS_PER_WARP;
619 const int num_warps = (num_rows + ROWS_PER_WARP - 1) / ROWS_PER_WARP;
620 const int num_blocks = (num_warps + WARPS_PER_TB - 1) / WARPS_PER_TB;
621
622 dim3 block_dim(WARP_SIZE_PARAM, WARPS_PER_TB);
623 topkGating<VPT, EXPERTS, WARPS_PER_TB, BYTES_PER_LDG, WARP_SIZE_PARAM, IndType, InputType, SF><<<num_blocks, block_dim, 0, stream>>>(
624 input, finished, output, num_rows, indices, source_row, k, start_expert, end_expert, renormalize, bias, routed_scaling_factor, is_padding);
625 }
626
627 #ifndef USE_ROCM
628 #define LAUNCH_TOPK(NUM_EXPERTS, WARPS_PER_TB, MAX_BYTES) \
629 static_assert(WARP_SIZE == 32, \
630 "Unsupported warp size. Only 32 is supported for CUDA"); \
631 topkGatingLauncherHelper<NUM_EXPERTS, WARPS_PER_TB, WARP_SIZE, MAX_BYTES, \
632 IndType, InputType, SF>( \
633 gating_output, nullptr, topk_weights, topk_indices, \
634 token_expert_indices, num_tokens, topk, 0, num_experts, renormalize, \
635 bias, routed_scaling_factor, stream, is_padding);
636 #else
637 #define LAUNCH_TOPK(NUM_EXPERTS, WARPS_PER_TB, MAX_BYTES) \
638 if (WARP_SIZE == 64) { \
639 topkGatingLauncherHelper<NUM_EXPERTS, WARPS_PER_TB, 64, MAX_BYTES, \
640 IndType, InputType, SF>( \
641 gating_output, nullptr, topk_weights, topk_indices, \
642 token_expert_indices, num_tokens, topk, 0, num_experts, renormalize, \
643 bias, routed_scaling_factor, stream, is_padding); \
644 } else if (WARP_SIZE == 32) { \
645 topkGatingLauncherHelper<NUM_EXPERTS, WARPS_PER_TB, 32, MAX_BYTES, \
646 IndType, InputType, SF>( \
647 gating_output, nullptr, topk_weights, topk_indices, \
648 token_expert_indices, num_tokens, topk, 0, num_experts, renormalize, \
649 bias, routed_scaling_factor, stream, is_padding); \
650 } else { \
651 assert(false && \
652 "Unsupported warp size. Only 32 and 64 are supported for ROCm"); \
653 }
654 #endif
655
656 template <typename IndType, typename InputType, ScoringFunc SF>
657 void topkGatingKernelLauncher(
658 const InputType* gating_output,
659 float* topk_weights,
660 IndType* topk_indices,
661 int* token_expert_indices,
662 float* workspace,
663 const int num_tokens,
664 const int num_experts,
665 const int topk,
666 const bool renormalize,
667 const float* bias,
668 const double routed_scaling_factor,
669 cudaStream_t stream,
670 const bool* is_padding) {
671 static constexpr int WARPS_PER_TB = 4;
672 static constexpr int BYTES_PER_LDG_POWER_OF_2 = 16;
673 #ifndef USE_ROCM
674 // for bfloat16 dtype, we need 4 bytes loading to make sure num_experts
675 // elements can be loaded by a warp
676 static constexpr int BYTES_PER_LDG_MULTIPLE_64 =
677 (std::is_same_v<InputType, __nv_bfloat16> || std::is_same_v<InputType, __half>) ? 4 : 8;
678 #endif
679 switch (num_experts) {
680 case 1:
681 LAUNCH_TOPK(1, WARPS_PER_TB, BYTES_PER_LDG_POWER_OF_2);
682 break;
683 case 2:
684 LAUNCH_TOPK(2, WARPS_PER_TB, BYTES_PER_LDG_POWER_OF_2);
685 break;
686 case 4:
687 LAUNCH_TOPK(4, WARPS_PER_TB, BYTES_PER_LDG_POWER_OF_2);
688 break;
689 case 8:
690 LAUNCH_TOPK(8, WARPS_PER_TB, BYTES_PER_LDG_POWER_OF_2);
691 break;
692 case 16:
693 LAUNCH_TOPK(16, WARPS_PER_TB, BYTES_PER_LDG_POWER_OF_2);
694 break;
695 case 32:
696 LAUNCH_TOPK(32, WARPS_PER_TB, BYTES_PER_LDG_POWER_OF_2);
697 break;
698 case 64:
699 LAUNCH_TOPK(64, WARPS_PER_TB, BYTES_PER_LDG_POWER_OF_2);
700 break;
701 case 128:
702 LAUNCH_TOPK(128, WARPS_PER_TB, BYTES_PER_LDG_POWER_OF_2);
703 break;
704 case 256:
705 LAUNCH_TOPK(256, WARPS_PER_TB, BYTES_PER_LDG_POWER_OF_2);
706 break;
707 case 512:
708 LAUNCH_TOPK(512, WARPS_PER_TB, BYTES_PER_LDG_POWER_OF_2);
709 break;
710 // (CUDA only) support multiples of 64 when num_experts is not power of 2.
711 // ROCm uses WARP_SIZE 64 so 8 bytes loading won't fit for some of num_experts,
712 // alternatively we can test 4 bytes loading and enable it in future.
713 #ifndef USE_ROCM
714 case 192:
715 LAUNCH_TOPK(192, WARPS_PER_TB, BYTES_PER_LDG_MULTIPLE_64);
716 break;
717 case 320:
718 LAUNCH_TOPK(320, WARPS_PER_TB, BYTES_PER_LDG_MULTIPLE_64);
719 break;
720 case 384:
721 LAUNCH_TOPK(384, WARPS_PER_TB, BYTES_PER_LDG_MULTIPLE_64);
722 break;
723 case 448:
724 LAUNCH_TOPK(448, WARPS_PER_TB, BYTES_PER_LDG_MULTIPLE_64);
725 break;
726 case 576:
727 LAUNCH_TOPK(576, WARPS_PER_TB, BYTES_PER_LDG_MULTIPLE_64);
728 break;
729 #endif
730 default: {
731 STD_TORCH_CHECK(workspace != nullptr,
732 "workspace must be provided for num_experts that are not a power of 2 or multiple of 64.");
733 static constexpr int TPB = 256;
734 if constexpr (SF == SCORING_SOFTMAX) {
735 moeSoftmax<TPB, InputType><<<num_tokens, TPB, 0, stream>>>(
736 gating_output, nullptr, workspace, num_experts);
737 } else if constexpr (SF == SCORING_SIGMOID) {
738 moeSigmoid<TPB, InputType><<<num_tokens, TPB, 0, stream>>>(
739 gating_output, nullptr, workspace, num_experts);
740 } else {
741 STD_TORCH_CHECK(false, "Unsupported scoring func");
742 }
743 moeTopK<TPB><<<num_tokens, TPB, 0, stream>>>(
744 workspace, nullptr, topk_weights, topk_indices, token_expert_indices,
745 num_experts, topk, 0, num_experts, renormalize, bias, routed_scaling_factor, is_padding);
746 }
747 }
748 }
749
750 } // namespace moe
751 } // namespace vllm
752
753
754 template<typename ComputeType, vllm::moe::ScoringFunc SF>
755 void dispatch_topk_launch(
756 torch::stable::Tensor& gating_output,
757 torch::stable::Tensor& topk_weights,
758 torch::stable::Tensor& topk_indices,
759 torch::stable::Tensor& token_expert_indices,
760 torch::stable::Tensor& softmax_workspace,
761 int num_tokens, int num_experts, int topk, bool renormalize,
762 std::optional<torch::stable::Tensor> bias,
763 double routed_scaling_factor,
764 cudaStream_t stream,
765 std::optional<torch::stable::Tensor> is_padding)
766 {
767 const float* bias_ptr = nullptr;
768 if (bias.has_value()) {
769 const torch::stable::Tensor& bias_tensor = bias.value();
770 STD_TORCH_CHECK(bias_tensor.scalar_type() == torch::headeronly::ScalarType::Float,
771 "bias tensor must be float32");
772 STD_TORCH_CHECK(bias_tensor.dim() == 1, "bias tensor must be 1D");
773 STD_TORCH_CHECK(bias_tensor.size(0) == num_experts,
774 "bias size mismatch, expected: ", num_experts);
775 STD_TORCH_CHECK(bias_tensor.is_contiguous(), "bias tensor must be contiguous");
776 bias_ptr = bias_tensor.const_data_ptr<float>();
777 }
778
779 const bool* is_padding_ptr = nullptr;
780 if (is_padding.has_value()) {
781 const torch::stable::Tensor& is_padding_tensor = is_padding.value();
782 STD_TORCH_CHECK(is_padding_tensor.scalar_type() == torch::headeronly::ScalarType::Bool,
783 "is_padding tensor must be bool");
784 STD_TORCH_CHECK(is_padding_tensor.dim() == 1, "is_padding tensor must be 1D");
785 STD_TORCH_CHECK(is_padding_tensor.size(0) == num_tokens,
786 "is_padding size mismatch, expected: ", num_tokens);
787 STD_TORCH_CHECK(is_padding_tensor.is_contiguous(), "is_padding tensor must be contiguous");
788 is_padding_ptr = is_padding_tensor.const_data_ptr<bool>();
789 }
790
791 if (topk_indices.scalar_type() == torch::headeronly::ScalarType::Int) {
792 vllm::moe::topkGatingKernelLauncher<int, ComputeType, SF>(
793 reinterpret_cast<const ComputeType*>(gating_output.const_data_ptr()),
794 topk_weights.mutable_data_ptr<float>(),
795 topk_indices.mutable_data_ptr<int>(),
796 token_expert_indices.mutable_data_ptr<int>(),
797 softmax_workspace.mutable_data_ptr<float>(),
798 num_tokens, num_experts, topk, renormalize,
799 bias_ptr, routed_scaling_factor, stream, is_padding_ptr);
800 } else if (topk_indices.scalar_type() == torch::headeronly::ScalarType::UInt32) {
801 vllm::moe::topkGatingKernelLauncher<uint32_t, ComputeType, SF>(
802 reinterpret_cast<const ComputeType*>(gating_output.const_data_ptr()),
803 topk_weights.mutable_data_ptr<float>(),
804 topk_indices.mutable_data_ptr<uint32_t>(),
805 token_expert_indices.mutable_data_ptr<int>(),
806 softmax_workspace.mutable_data_ptr<float>(),
807 num_tokens, num_experts, topk, renormalize,
808 bias_ptr, routed_scaling_factor, stream, is_padding_ptr);
809 } else {
810 STD_TORCH_CHECK(topk_indices.scalar_type() == torch::headeronly::ScalarType::Long);
811 vllm::moe::topkGatingKernelLauncher<int64_t, ComputeType, SF>(
812 reinterpret_cast<const ComputeType*>(gating_output.const_data_ptr()),
813 topk_weights.mutable_data_ptr<float>(),
814 topk_indices.mutable_data_ptr<int64_t>(),
815 token_expert_indices.mutable_data_ptr<int>(),
816 softmax_workspace.mutable_data_ptr<float>(),
817 num_tokens, num_experts, topk, renormalize,
818 bias_ptr, routed_scaling_factor, stream, is_padding_ptr);
819 }
820 }
821
822 void topk_softmax(
823 torch::stable::Tensor& topk_weights, // [num_tokens, topk]
824 torch::stable::Tensor& topk_indices, // [num_tokens, topk]
825 torch::stable::Tensor& token_expert_indices, // [num_tokens, topk]
826 torch::stable::Tensor& gating_output, // [num_tokens, num_experts]
827 bool renormalize,
828 std::optional<torch::stable::Tensor> bias,
829 std::optional<torch::stable::Tensor> is_padding)
830 {
831 const int num_experts = gating_output.size(-1);
832 const auto num_tokens = gating_output.numel() / num_experts;
833 const int topk = topk_weights.size(-1);
834
835 const bool is_pow_2 = (num_experts != 0) && ((num_experts & (num_experts - 1)) == 0);
836 const bool needs_workspace = !is_pow_2 || num_experts > 256;
837 const int64_t workspace_size = needs_workspace ? num_tokens * num_experts : 0;
838
839 torch::stable::accelerator::DeviceGuard guard(gating_output.get_device_index());
840 const cudaStream_t stream =
841 get_current_cuda_stream(gating_output.get_device_index());
842 auto softmax_workspace = torch::stable::new_empty(
843 gating_output, {workspace_size}, torch::headeronly::ScalarType::Float);
844
845 if (gating_output.scalar_type() == torch::headeronly::ScalarType::Float) {
846 dispatch_topk_launch<float, vllm::moe::SCORING_SOFTMAX>(gating_output, topk_weights, topk_indices,
847 token_expert_indices, softmax_workspace, num_tokens, num_experts, topk, renormalize,
848 bias, 1.0, stream, is_padding);
849 } else if (gating_output.scalar_type() == torch::headeronly::ScalarType::Half) {
850 dispatch_topk_launch<__half, vllm::moe::SCORING_SOFTMAX>(gating_output, topk_weights, topk_indices,
851 token_expert_indices, softmax_workspace, num_tokens, num_experts, topk, renormalize,
852 bias, 1.0, stream, is_padding);
853 } else if (gating_output.scalar_type() == torch::headeronly::ScalarType::BFloat16) {
854 dispatch_topk_launch<__nv_bfloat16, vllm::moe::SCORING_SOFTMAX>(gating_output, topk_weights, topk_indices,
855 token_expert_indices, softmax_workspace, num_tokens, num_experts, topk, renormalize,
856 bias, 1.0, stream, is_padding);
857 } else {
858 STD_TORCH_CHECK(false, "Unsupported gating_output data type: ", gating_output.scalar_type());
859 }
860 }
861
862 void topk_sigmoid(
863 torch::stable::Tensor& topk_weights, // [num_tokens, topk]
864 torch::stable::Tensor& topk_indices, // [num_tokens, topk]
865 torch::stable::Tensor& token_expert_indices, // [num_tokens, topk]
866 torch::stable::Tensor& gating_output, // [num_tokens, num_experts]
867 bool renormalize,
868 std::optional<torch::stable::Tensor> bias,
869 double routed_scaling_factor,
870 std::optional<torch::stable::Tensor> is_padding)
871 {
872 const int num_experts = gating_output.size(-1);
873 const auto num_tokens = gating_output.numel() / num_experts;
874 const int topk = topk_weights.size(-1);
875
876 const bool is_pow_2 = (num_experts != 0) && ((num_experts & (num_experts - 1)) == 0);
877 const bool needs_workspace = !is_pow_2 || num_experts > 256;
878 const int64_t workspace_size = needs_workspace ? num_tokens * num_experts : 0;
879
880 torch::stable::accelerator::DeviceGuard guard(gating_output.get_device_index());
881 const cudaStream_t stream =
882 get_current_cuda_stream(gating_output.get_device_index());
883 auto workspace = torch::stable::new_empty(
884 gating_output, {workspace_size}, torch::headeronly::ScalarType::Float);
885
886 if (gating_output.scalar_type() == torch::headeronly::ScalarType::Float) {
887 dispatch_topk_launch<float, vllm::moe::SCORING_SIGMOID>(gating_output, topk_weights, topk_indices,
888 token_expert_indices, workspace, num_tokens, num_experts, topk, renormalize,
889 bias, routed_scaling_factor, stream, is_padding);
890 } else if (gating_output.scalar_type() == torch::headeronly::ScalarType::Half) {
891 dispatch_topk_launch<__half, vllm::moe::SCORING_SIGMOID>(gating_output, topk_weights, topk_indices,
892 token_expert_indices, workspace, num_tokens, num_experts, topk, renormalize,
893 bias, routed_scaling_factor, stream, is_padding);
894 } else if (gating_output.scalar_type() == torch::headeronly::ScalarType::BFloat16) {
895 dispatch_topk_launch<__nv_bfloat16, vllm::moe::SCORING_SIGMOID>(gating_output, topk_weights, topk_indices,
896 token_expert_indices, workspace, num_tokens, num_experts, topk, renormalize,
897 bias, routed_scaling_factor, stream, is_padding);
898 } else {
899 STD_TORCH_CHECK(false, "Unsupported gating_output data type: ", gating_output.scalar_type());
900 }
901 }