main
cu 865 lines 31.3 KB
Raw
1 /*
2 * Adapted from
3 * https://github.com/sgl-project/sglang/blob/main/sgl-kernel/csrc/gemm/dsv3_fused_a_gemm.cu
4 * which was adapted from
5 * https://github.com/NVIDIA/TensorRT-LLM/blob/619709fc33bd5dc268f19d6a741fe7ed51c0f8f5/cpp/tensorrt_llm/kernels/dsv3MinLatencyKernels/dsv3FusedAGemm.cu
6 *
7 * Copyright (c) 2019-2024, NVIDIA CORPORATION. All rights reserved.
8 * Copyright (c) 2021, NAVER Corp. Authored by CLOVA.
9 *
10 * Licensed under the Apache License, Version 2.0 (the "License");
11 * you may not use this file except in compliance with the License.
12 * You may obtain a copy of the License at
13 *
14 * http://www.apache.org/licenses/LICENSE-2.0
15 *
16 * Unless required by applicable law or agreed to in writing, software
17 * distributed under the License is distributed on an "AS IS" BASIS,
18 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
19 * See the License for the specific language governing permissions and
20 * limitations under the License.
21 */
22
23 #include <torch/csrc/stable/library.h>
24 #include <torch/csrc/stable/tensor.h>
25 #include <torch/headeronly/core/ScalarType.h>
26
27 #include "core/registration.h"
28 #include "libtorch_stable/torch_utils.h"
29
30 #include <cuda_bf16.h>
31 #include <cuda_runtime.h>
32
33 #include <cstdlib>
34 #include <mutex>
35
36 namespace {
37
38 inline int getSMVersion() {
39 auto* props = get_device_prop();
40 return props->major * 10 + props->minor;
41 }
42
43 inline bool getEnvEnablePDL() {
44 static std::once_flag flag;
45 static bool enablePDL = false;
46 std::call_once(flag, [&]() {
47 if (getSMVersion() >= 90) {
48 char const* env = std::getenv("TRTLLM_ENABLE_PDL");
49 enablePDL = env && env[0] == '1' && env[1] == '\0';
50 }
51 });
52 return enablePDL;
53 }
54
55 } // namespace
56
57 using bf16_t = __nv_bfloat16;
58
59 __device__ void hmma_16_8_16_f32acc_bf16ab(float (&d_reg)[4],
60 const bf16_t (&a_reg)[8],
61 const bf16_t (&b_reg)[4],
62 float const (&c_reg)[4]) {
63 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
64 uint32_t a0 = *reinterpret_cast<uint32_t const*>(a_reg + 0);
65 uint32_t a1 = *reinterpret_cast<uint32_t const*>(a_reg + 2);
66 uint32_t a2 = *reinterpret_cast<uint32_t const*>(a_reg + 4);
67 uint32_t a3 = *reinterpret_cast<uint32_t const*>(a_reg + 6);
68 uint32_t b0 = *reinterpret_cast<uint32_t const*>(b_reg + 0);
69 uint32_t b1 = *reinterpret_cast<uint32_t const*>(b_reg + 2);
70 asm volatile(
71 "mma.sync.aligned.m16n8k16.row.col.f32.bf16.bf16.f32 "
72 "{%0, %1, %2, %3},"
73 "{%4, %5, %6, %7},"
74 "{%8, %9},"
75 "{%10, %11, %12, %13};\n"
76 : "=f"(d_reg[0]), "=f"(d_reg[1]), "=f"(d_reg[2]), "=f"(d_reg[3])
77 : "r"(a0), "r"(a1), "r"(a2), "r"(a3), "r"(b0), "r"(b1), "f"(d_reg[0]),
78 "f"(d_reg[1]), "f"(d_reg[2]), "f"(d_reg[3]));
79 #endif
80 }
81
82 extern "C" {
83 __device__ uint32_t __nvvm_get_smem_pointer(void*);
84 }
85
86 __device__ void ldgsts_128(void const* gPtr, void* sPtr, uint32_t pred) {
87 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
88 if (pred) {
89 uint32_t smemPtrAsUint32 = __nvvm_get_smem_pointer(sPtr);
90 asm volatile("cp.async.cg.shared.global.L2::128B [%0], [%1], %2;\n" ::"r"(
91 smemPtrAsUint32),
92 "l"(gPtr), "n"(16));
93 }
94 #endif
95 }
96
97 __device__ void ldsm_x4(void* smem_ptr, uint32_t* reg_ptr) {
98 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
99 asm volatile(
100 "ldmatrix.sync.aligned.x4.m8n8.shared.b16 {%0, %1, %2, %3}, [%4];\n"
101 : "=r"(reg_ptr[0]), "=r"(reg_ptr[1]), "=r"(reg_ptr[2]), "=r"(reg_ptr[3])
102 : "r"(__nvvm_get_smem_pointer(smem_ptr)));
103 #endif
104 }
105
106 template <class Type>
107 __device__ int apply_swizzle_343_on_elem_row_col(int row_idx_, int col_idx_) {
108 uint32_t row_idx = *reinterpret_cast<uint32_t*>(&row_idx_);
109 uint32_t col_idx = *reinterpret_cast<uint32_t*>(&col_idx_);
110 row_idx = row_idx % 8;
111 row_idx = row_idx * (16 / sizeof(Type));
112 col_idx = col_idx ^ row_idx;
113 return *reinterpret_cast<int*>(&col_idx);
114 }
115
116 __device__ void initialize_barrier(
117 uint64_t* smem_barrier, // 64 bits user-manged barrier in smem
118 int thread_count =
119 1) // Thread count expected to arrive/wait on this barrier
120 {
121 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
122 uint32_t smem_int_ptr = __nvvm_get_smem_pointer(smem_barrier);
123 asm volatile("mbarrier.init.shared::cta.b64 [%0], %1;\n" ::"r"(smem_int_ptr),
124 "r"(thread_count));
125 #endif
126 }
127
128 // Barrier wait
129 __device__ void wait_barrier(
130 uint64_t* smem_barrier, // 64 bits user-manged barrier in smem
131 int phase_bit) // Current phase bit the barrier waiting to flip
132 {
133 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
134 uint32_t smem_int_ptr = __nvvm_get_smem_pointer(smem_barrier);
135 asm volatile(
136 "{\n"
137 ".reg .pred P1;\n"
138 "LAB_WAIT:\n"
139 "mbarrier.try_wait.parity.shared::cta.b64 P1, [%0], %1;\n"
140 "@P1 bra DONE;\n"
141 "bra LAB_WAIT;\n"
142 "DONE:\n"
143 "}\n" ::"r"(smem_int_ptr),
144 "r"(phase_bit));
145 #endif
146 }
147
148 __device__ bool try_wait_barrier(uint64_t* smem_ptr, int phase_bit) {
149 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
150 uint32_t wait_complete;
151 uint32_t smem_int_ptr = __nvvm_get_smem_pointer(smem_ptr);
152 asm volatile(
153 "{\n\t"
154 ".reg .pred P1; \n\t"
155 "mbarrier.try_wait.parity.shared::cta.b64 P1, [%1], %2; \n\t"
156 "selp.b32 %0, 1, 0, P1; \n\t"
157 "}"
158 : "=r"(wait_complete)
159 : "r"(smem_int_ptr), "r"(phase_bit));
160 return static_cast<bool>(wait_complete);
161 #endif
162 return false;
163 }
164
165 // Barrier arrive
166 __device__ void arrive_barrier(
167 uint64_t* smem_barrier) // 64 bits user-manged barrier in smem
168 {
169 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
170 uint32_t smem_int_ptr = __nvvm_get_smem_pointer(smem_barrier);
171 asm volatile(
172 "{\n"
173 ".reg .b64 state; \n"
174 "mbarrier.arrive.shared::cta.b64 state, [%0];\n"
175 "}\n" ::"r"(smem_int_ptr));
176 #endif
177 }
178
179 __device__ void ldgsts_arrive(uint64_t* smem_barrier) {
180 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
181 uint32_t smem_int_ptr = __nvvm_get_smem_pointer(smem_barrier);
182 asm volatile("cp.async.mbarrier.arrive.noinc.shared.b64 [%0];"
183 :
184 : "r"(smem_int_ptr));
185 #endif
186 }
187
188 template <int gemm_k, int tile_m, int tile_k, int stage_cnt>
189 struct GmemLoaderA {
190 static constexpr int elem_bytes = 2;
191 static constexpr int vec_bytes = 16;
192 static constexpr int vec_elems = vec_bytes / elem_bytes;
193 static constexpr int thread_cnt = 64;
194 static_assert((tile_m * tile_k) % (vec_elems * thread_cnt) == 0);
195 static constexpr int a_inst_cnt_per_iter =
196 (tile_m * tile_k) / (vec_elems * thread_cnt);
197 static_assert(gemm_k % tile_k == 0);
198 static constexpr int k_iter_cnt = gemm_k / tile_k;
199
200 // Extra params to keep the order of k reduction...
201 static constexpr int mma_warp_cnt = 4;
202 static constexpr int per_mma_warp_k = tile_k / mma_warp_cnt;
203 static constexpr int k_each_chunk = gemm_k / mma_warp_cnt;
204
205 private:
206 __device__ int k_project(int tile_k_idx) {
207 return (tile_k_idx / per_mma_warp_k * k_each_chunk) +
208 (tile_k_idx % per_mma_warp_k);
209 }
210
211 public:
212 __device__ GmemLoaderA(bf16_t const* gmem_a_local_, bf16_t* smem_a_,
213 uint64_t* smem_barrier_)
214 : gmem_a(gmem_a_local_),
215 smem_a(smem_a_),
216 smem_barrier(smem_barrier_),
217 local_tid(threadIdx.x % thread_cnt) {}
218
219 __device__ void prepare() {
220 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
221 // swizzle, that's what we want.
222 #pragma unroll
223 for (int i = 0; i < a_inst_cnt_per_iter; i++) {
224 int linear_idx = local_tid * vec_elems + i * thread_cnt * vec_elems;
225 int m_idx = linear_idx / tile_k;
226 int k_idx = linear_idx % tile_k;
227 k_idx = apply_swizzle_343_on_elem_row_col<bf16_t>(m_idx, k_idx);
228 a_smem_offsets[i] = m_idx * tile_k + k_idx;
229 }
230 #endif
231 }
232
233 __device__ void issue_mainloop() {
234 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
235 #pragma unroll 1
236 for (int loop_idx = 0; loop_idx < k_iter_cnt; loop_idx++) {
237 if (need_wait) {
238 wait_barrier(smem_barrier + 1 + stage_idx * 2, phase_bit);
239 }
240 int next_stage_idx = stage_idx + 1;
241 int next_phase_bit =
242 next_stage_idx == stage_cnt ? phase_bit ^ 1 : phase_bit;
243 next_stage_idx = next_stage_idx == stage_cnt ? 0 : next_stage_idx;
244 if (loop_idx != k_iter_cnt - 1) {
245 need_wait = !try_wait_barrier(smem_barrier + 1 + next_stage_idx * 2,
246 next_phase_bit);
247 }
248
249 #pragma unroll
250 for (int i = 0; i < a_inst_cnt_per_iter; i++) {
251 int smem_offset = a_smem_offsets[i];
252 bf16_t* smem_ptr_this_iter =
253 smem_a + stage_idx * tile_m * tile_k + smem_offset;
254 int linear_idx = local_tid * vec_elems + i * thread_cnt * vec_elems;
255 int m_idx = linear_idx / tile_k;
256 int k_idx = linear_idx % tile_k;
257 int gmem_offset = m_idx * gemm_k + k_project(k_idx);
258 bf16_t const* gmem_ptr_this_iter = gmem_a + gmem_offset;
259 ldgsts_128(gmem_ptr_this_iter, smem_ptr_this_iter, true);
260 }
261 ldgsts_arrive(smem_barrier + stage_idx * 2);
262
263 stage_idx = next_stage_idx;
264 phase_bit = next_phase_bit;
265 gmem_a += per_mma_warp_k;
266 }
267 #endif
268 }
269
270 bf16_t const* gmem_a;
271 bf16_t* smem_a;
272 uint64_t* smem_barrier;
273 int local_tid;
274 int stage_idx = 0;
275 int phase_bit = 1;
276 bool need_wait = true;
277
278 // per smem_stage, store with swizzle information
279 int a_smem_offsets[a_inst_cnt_per_iter];
280 };
281
282 template <int gemm_k, int tile_n, int tile_k, int stage_cnt>
283 struct GmemLoaderB {
284 static constexpr int elem_bytes = 2;
285 static constexpr int vec_bytes = 16;
286 static constexpr int vec_elems = vec_bytes / elem_bytes;
287 static constexpr int thread_cnt = 64;
288 static_assert((tile_n * tile_k) % (vec_elems * thread_cnt) == 0);
289 static constexpr int b_inst_cnt_per_iter =
290 (tile_n * tile_k) / (vec_elems * thread_cnt);
291 static_assert(gemm_k % tile_k == 0);
292 static constexpr int k_iter_cnt = gemm_k / tile_k;
293
294 // Extra params to keep the order of k reduction...
295 static constexpr int mma_warp_cnt = 4;
296 static constexpr int per_mma_warp_k = tile_k / mma_warp_cnt;
297 static constexpr int k_each_chunk = gemm_k / mma_warp_cnt;
298
299 private:
300 __device__ int k_project(int tile_k_idx) {
301 return (tile_k_idx / per_mma_warp_k * k_each_chunk) +
302 (tile_k_idx % per_mma_warp_k);
303 }
304
305 public:
306 __device__ GmemLoaderB(bf16_t const* gmem_b_local_, bf16_t* smem_b_,
307 uint64_t* smem_barrier_, int gemm_n_)
308 : gmem_b(gmem_b_local_),
309 smem_b(smem_b_),
310 smem_barrier(smem_barrier_),
311 gemm_n(gemm_n_),
312 local_tid(threadIdx.x % thread_cnt) {}
313
314 __device__ void prepare() {
315 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
316 // swizzle, that's what we want.
317 #pragma unroll
318 for (int i = 0; i < b_inst_cnt_per_iter; i++) {
319 int linear_idx = local_tid * vec_elems + i * thread_cnt * vec_elems;
320 int n_idx = linear_idx / tile_k;
321 int k_idx = linear_idx % tile_k;
322 k_idx = apply_swizzle_343_on_elem_row_col<bf16_t>(n_idx, k_idx);
323 b_smem_offsets[i] = n_idx * tile_k + k_idx;
324 preds[i] = n_idx < gemm_n;
325 }
326 #endif
327 }
328
329 __device__ void issue_mainloop(int row_stride) {
330 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
331 cudaGridDependencySynchronize();
332 #pragma unroll 1
333 for (int loop_idx = 0; loop_idx < k_iter_cnt; loop_idx++) {
334 if (need_wait) {
335 wait_barrier(smem_barrier + 1 + stage_idx * 2, phase_bit);
336 }
337 int next_stage_idx = stage_idx + 1;
338 int next_phase_bit =
339 next_stage_idx == stage_cnt ? phase_bit ^ 1 : phase_bit;
340 next_stage_idx = next_stage_idx == stage_cnt ? 0 : next_stage_idx;
341 if (loop_idx != k_iter_cnt - 1) {
342 need_wait = !try_wait_barrier(smem_barrier + 1 + next_stage_idx * 2,
343 next_phase_bit);
344 }
345 #pragma unroll
346 for (int i = 0; i < b_inst_cnt_per_iter; i++) {
347 int smem_offset = b_smem_offsets[i];
348 bf16_t* smem_ptr_this_iter =
349 smem_b + stage_idx * tile_n * tile_k + smem_offset;
350 int linear_idx = local_tid * vec_elems + i * thread_cnt * vec_elems;
351 int n_idx = linear_idx / tile_k;
352 int k_idx = linear_idx % tile_k;
353 int gmem_offset = n_idx * row_stride + k_project(k_idx);
354 bf16_t const* gmem_ptr_this_iter = gmem_b + gmem_offset;
355 ldgsts_128(gmem_ptr_this_iter, smem_ptr_this_iter, preds[i]);
356 }
357 ldgsts_arrive(smem_barrier + stage_idx * 2);
358
359 stage_idx = next_stage_idx;
360 phase_bit = next_phase_bit;
361 gmem_b += per_mma_warp_k;
362 }
363 #endif
364 }
365
366 bf16_t const* gmem_b;
367 bf16_t* smem_b;
368 uint64_t* smem_barrier;
369 int gemm_n;
370 int local_tid;
371 int stage_idx = 0;
372 int phase_bit = 1;
373 bool need_wait = true;
374
375 // per smem_stage, store with swizzle information
376 int b_smem_offsets[b_inst_cnt_per_iter];
377 uint32_t preds[b_inst_cnt_per_iter];
378 };
379
380 template <int gemm_m, int gemm_k, int tile_m, int tile_n, int tile_k,
381 int stage_cnt>
382 struct MmaComputer {
383 static constexpr int elem_bytes = 2;
384 static constexpr int thread_cnt = 128;
385 static_assert(gemm_k % tile_k == 0);
386 static_assert(tile_k % (thread_cnt / 32) == 0);
387 static constexpr int per_warp_tile_k = tile_k / (thread_cnt / 32);
388 static constexpr int k_iter_cnt = gemm_k / tile_k;
389 static constexpr int k_phase_cnt = per_warp_tile_k / 16;
390 static constexpr int m_iter_cnt = (tile_m + 15) / 16;
391 static constexpr int n_iter_cnt =
392 (tile_n + 7) /
393 8; // Possible to have non-1 n_iter_cnt for ab_swap m16 case.
394 static_assert(m_iter_cnt == 1 || m_iter_cnt == 2);
395 static_assert(n_iter_cnt == 1 || n_iter_cnt == 2);
396
397 __device__ MmaComputer(bf16_t* gmem_c_local_, bf16_t* smem_a_,
398 bf16_t* smem_b_, uint64_t* smem_barrier_,
399 int warp_idx_, int gemm_n_)
400 : gmem_c(gmem_c_local_),
401 smem_a(smem_a_),
402 smem_b(smem_b_),
403 smem_barrier(smem_barrier_),
404 warp_idx(warp_idx_ - (thread_cnt / 32)),
405 gemm_n(gemm_n_) {}
406
407 private:
408 __device__ constexpr int internal_b_atom_func(int tid) {
409 if constexpr (tile_n < 8) {
410 return (tid % tile_n) + ((tid % 8) / tile_n * 0) + tid / 8 * 8 * tile_n;
411 } else {
412 return (tid % 8) + ((tid % 32) / 8 * (tile_n * 8));
413 }
414 }
415
416 public:
417 __device__ void prepare() {
418 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
419 #pragma unroll
420 for (int m = 0; m < m_iter_cnt; m++) {
421 #pragma unroll
422 for (int i = 0; i < k_phase_cnt; i++) {
423 int linear_idx = (lane_idx % 16) + (lane_idx / 16) * 128 + i * 256;
424 int m_idx = linear_idx % 16 + m * 16;
425 int k_idx = linear_idx / 16 + warp_k_offset_in_tile_k;
426 k_idx = apply_swizzle_343_on_elem_row_col<bf16_t>(m_idx, k_idx);
427 a_smem_offsets[m][i] = m_idx * tile_k + k_idx;
428 }
429 }
430 #pragma unroll
431 for (int n_iter_idx = 0; n_iter_idx < n_iter_cnt; n_iter_idx++) {
432 #pragma unroll
433 for (int i = 0; i < k_phase_cnt; i += 2) { // Special i+=2 for B.
434 int linear_idx =
435 internal_b_atom_func(lane_idx) + i * tile_n * 16 + n_iter_idx * 8;
436 int n_idx = linear_idx % tile_n;
437 int k_idx = linear_idx / tile_n + warp_k_offset_in_tile_k;
438 k_idx = apply_swizzle_343_on_elem_row_col<bf16_t>(n_idx, k_idx);
439 b_smem_offsets[n_iter_idx][i] = n_idx * tile_k + k_idx;
440 }
441 }
442 #endif
443 }
444
445 __device__ void issue_mainloop() {
446 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
447 #pragma unroll 1
448 for (int loop_idx = 0; loop_idx < k_iter_cnt; loop_idx++) {
449 wait_barrier(smem_barrier + 0 + stage_idx * 2, phase_bit);
450
451 #pragma unroll
452 for (int m = 0; m < m_iter_cnt; m++) {
453 #pragma unroll
454 for (int i = 0; i < k_phase_cnt; i++) {
455 int smem_offset = a_smem_offsets[m][i];
456 bf16_t* smem_ptr_this_iter =
457 smem_a + stage_idx * tile_m * tile_k + smem_offset;
458 ldsm_x4(smem_ptr_this_iter, reinterpret_cast<uint32_t*>(a_reg[m][i]));
459 }
460 }
461
462 #pragma unroll
463 for (int n_iter_idx = 0; n_iter_idx < n_iter_cnt; n_iter_idx++) {
464 #pragma unroll
465 for (int i = 0; i < k_phase_cnt; i += 2) {
466 int smem_offset = b_smem_offsets[n_iter_idx][i];
467 bf16_t* smem_ptr_this_iter =
468 smem_b + stage_idx * tile_n * tile_k + smem_offset;
469 ldsm_x4(smem_ptr_this_iter,
470 reinterpret_cast<uint32_t*>(b_reg[n_iter_idx][i]));
471 }
472 }
473
474 #pragma unroll
475 for (int k_iter_idx = 0; k_iter_idx < k_phase_cnt; k_iter_idx++) {
476 #pragma unroll
477 for (int n_iter_idx = 0; n_iter_idx < n_iter_cnt; n_iter_idx++) {
478 #pragma unroll
479 for (int m = 0; m < m_iter_cnt; m++) {
480 hmma_16_8_16_f32acc_bf16ab(
481 acc_reg[m][n_iter_idx], a_reg[m][k_iter_idx],
482 b_reg[n_iter_idx][k_iter_idx], acc_reg[m][n_iter_idx]);
483 }
484 }
485 }
486 ::arrive_barrier(smem_barrier + 1 + stage_idx * 2);
487 stage_idx += 1;
488 phase_bit = stage_idx == stage_cnt ? phase_bit ^ 1 : phase_bit;
489 stage_idx = stage_idx == stage_cnt ? 0 : stage_idx;
490 }
491 #endif
492 }
493
494 __device__ void epi() {
495 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
496 asm volatile("bar.sync %0, %1;" : : "r"(1), "r"(thread_cnt));
497 // reorganize the acc_reg
498 constexpr int thread_m = 2 * m_iter_cnt;
499 constexpr int thread_n = 2 * n_iter_cnt;
500 constexpr int cta_mma_n = n_iter_cnt * 8;
501 float acc_reg_reorg[thread_m][thread_n];
502
503 for (int i = 0; i < thread_m; i++) {
504 for (int j = 0; j < thread_n; j++) {
505 acc_reg_reorg[i][j] = acc_reg[i / 2][j / 2][(j % 2) + (i % 2) * 2];
506 }
507 }
508
509 // 4 x cosize(smem_c_layout)
510 float* smem_c = reinterpret_cast<float*>(smem_a);
511 // coord -> index
512 auto smem_c_index_func = [&](int m_idx, int n_idx) {
513 int group_rows = 32 / cta_mma_n;
514 int group_cnt = 2;
515 return (m_idx % group_rows * cta_mma_n) +
516 (m_idx / group_rows * (32 + group_cnt)) + n_idx;
517 };
518 constexpr int cosize_smem_c = ((tile_m * cta_mma_n) / 32) * (32 + 2);
519
520 // This should be optimized to STS.64 but can not be STS.128 due to the bank
521 // index.
522 #pragma unroll
523 for (int m_idx_thread = 0; m_idx_thread < thread_m; m_idx_thread++) {
524 #pragma unroll
525 for (int n_idx_thread = 0; n_idx_thread < thread_n; n_idx_thread++) {
526 int m_idx =
527 (lane_idx / 4) + (m_idx_thread % 2) * 8 + (m_idx_thread / 2) * 16;
528 int n_idx =
529 ((lane_idx % 4) * 2) + (n_idx_thread % 2) + (n_idx_thread / 2) * 8;
530 smem_c[cosize_smem_c * warp_idx + smem_c_index_func(m_idx, n_idx)] =
531 acc_reg_reorg[m_idx_thread][n_idx_thread];
532 }
533 }
534 asm volatile("bar.sync %0, %1;" : : "r"(1), "r"(thread_cnt));
535
536 if (warp_idx == 0) {
537 constexpr int final_acc_reg_cnt = (tile_m * tile_n + 31) / 32;
538 float acc_final[final_acc_reg_cnt]{};
539
540 #pragma unroll
541 for (int reg_idx = 0; reg_idx < final_acc_reg_cnt; reg_idx++) {
542 int linear_idx = reg_idx * 32 + lane_idx;
543 int m_idx = linear_idx % tile_m;
544 int n_idx = linear_idx / tile_m;
545 acc_final[reg_idx] +=
546 smem_c[smem_c_index_func(m_idx, n_idx) + 0 * cosize_smem_c] +
547 smem_c[smem_c_index_func(m_idx, n_idx) + 1 * cosize_smem_c] +
548 smem_c[smem_c_index_func(m_idx, n_idx) + 2 * cosize_smem_c] +
549 smem_c[smem_c_index_func(m_idx, n_idx) + 3 * cosize_smem_c];
550 }
551
552 #pragma unroll
553 for (int reg_idx = 0; reg_idx < final_acc_reg_cnt; reg_idx++) {
554 int linear_idx = reg_idx * 32 + lane_idx;
555 int m_idx = linear_idx % tile_m;
556 int n_idx = linear_idx / tile_m;
557 if (m_idx < tile_m && n_idx < gemm_n) {
558 gmem_c[n_idx * gemm_m + m_idx] = acc_final[reg_idx];
559 }
560 }
561 }
562 #endif
563 }
564
565 bf16_t* gmem_c;
566 bf16_t* smem_a;
567 bf16_t* smem_b;
568 uint64_t* smem_barrier;
569 int warp_idx;
570 int gemm_n;
571 int stage_idx = 0;
572 int phase_bit = 0;
573 int lane_idx = threadIdx.x % 32;
574 int warp_k_offset_in_tile_k = warp_idx * per_warp_tile_k;
575
576 int a_smem_offsets[m_iter_cnt][k_phase_cnt];
577 int b_smem_offsets[n_iter_cnt][k_phase_cnt];
578
579 bf16_t a_reg[m_iter_cnt][k_phase_cnt][8];
580 bf16_t b_reg[n_iter_cnt][k_phase_cnt][4];
581 float acc_reg[m_iter_cnt][n_iter_cnt][4]{};
582 };
583
584 // AB swapped, kernel is k-major, k-major, m-major
585 template <int batch_size, int gemm_m, int gemm_k, int tile_m, int tile_n,
586 int tile_k, int stage_cnt>
587 __global__ __launch_bounds__(256, 1) void fused_a_gemm_kernel(
588 bf16_t* output, bf16_t const* mat_a, bf16_t const* mat_b, int gemm_n,
589 int input_row_stride) {
590 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
591 constexpr int load_thread_cnt = 128;
592 constexpr int compute_thread_cnt = 128;
593 constexpr int thread_cnt = load_thread_cnt + compute_thread_cnt;
594 (void)thread_cnt;
595 static_assert(gemm_m % 16 == 0);
596 static_assert(gemm_k % tile_k == 0);
597 static_assert(gemm_m % tile_m == 0);
598 static_assert(
599 tile_k == 128 || tile_k == 256 || tile_k == 512 ||
600 tile_k == 1024); // tile_k must be larger than 64 since 4 warp splitK.
601 static_assert(tile_m == 16 || tile_m == 32);
602 constexpr int g2s_vec_bytes = 16;
603 constexpr int a_elem_bytes = 2;
604 constexpr int b_elem_bytes = 2;
605 static_assert((tile_m * a_elem_bytes + tile_n * b_elem_bytes) * tile_k *
606 stage_cnt <=
607 225 * 1024);
608 static_assert((tile_m * tile_k * a_elem_bytes) %
609 (load_thread_cnt * g2s_vec_bytes) ==
610 0);
611 static_assert((tile_n * tile_k * b_elem_bytes) %
612 (load_thread_cnt * g2s_vec_bytes) ==
613 0);
614
615 extern __shared__ char smem[];
616 uint64_t* smem_barrier = reinterpret_cast<uint64_t*>(
617 smem); // producer,consumer; producer,consumer; ...
618 bf16_t* smem_a = reinterpret_cast<bf16_t*>(smem + (stage_cnt * 8 * 2 + 1024) /
619 1024 * 1024);
620 bf16_t* smem_b = smem_a + tile_m * tile_k * stage_cnt;
621
622 int cta_m_idx = tile_m * blockIdx.x;
623 int cta_n_idx = tile_n * blockIdx.y;
624 bf16_t const* gmem_a_local = mat_a + cta_m_idx * gemm_k;
625 bf16_t const* gmem_b_local = mat_b + cta_n_idx * input_row_stride;
626 bf16_t* gmem_c_local = output + cta_n_idx * gemm_m + cta_m_idx;
627
628 int warp_idx = __shfl_sync(0xffffffff, threadIdx.x / 32, 0);
629
630 if (warp_idx == 4) {
631 for (int i = 0; i < stage_cnt; i++) {
632 initialize_barrier(smem_barrier + i * 2 + 0,
633 load_thread_cnt); // producer
634 initialize_barrier(smem_barrier + i * 2 + 1,
635 compute_thread_cnt); // consumer
636 }
637 }
638 __syncthreads();
639
640 if (warp_idx < 2) {
641 GmemLoaderA<gemm_k, tile_m, tile_k, stage_cnt> a_loader(
642 gmem_a_local, smem_a, smem_barrier);
643 a_loader.prepare();
644 a_loader.issue_mainloop();
645 } else if (warp_idx < 4) {
646 GmemLoaderB<gemm_k, tile_n, tile_k, stage_cnt> b_loader(
647 gmem_b_local, smem_b, smem_barrier, gemm_n);
648 b_loader.prepare();
649 b_loader.issue_mainloop(input_row_stride);
650 } else {
651 MmaComputer<gemm_m, gemm_k, tile_m, tile_n, tile_k, stage_cnt> mma_computer(
652 gmem_c_local, smem_a, smem_b, smem_barrier, warp_idx, gemm_n);
653 mma_computer.prepare();
654 mma_computer.issue_mainloop();
655 mma_computer.epi();
656 }
657 cudaTriggerProgrammaticLaunchCompletion();
658 #endif
659 }
660
661 template <typename T, int kHdIn, int kHdOut, int kTileN, int kTileK = 256,
662 int kTileM = 16>
663 void invokeFusedAGemm(T* output, T const* mat_a, T const* mat_b, int num_tokens,
664 int input_row_stride, cudaStream_t const stream,
665 bool enable_pdl) {
666 constexpr int gemm_m = kHdOut;
667 int const gemm_n = num_tokens;
668 constexpr int gemm_k = kHdIn;
669 constexpr int batch_size = 1;
670 std::swap(mat_a, mat_b);
671 constexpr int tile_m = kTileM;
672 constexpr int tile_n = kTileN;
673 constexpr int tile_k = kTileK;
674 constexpr int max_stage_cnt =
675 1024 * 192 / ((tile_m + tile_n) * tile_k * sizeof(bf16_t));
676 constexpr int k_iter_cnt = gemm_k / tile_k;
677 constexpr int stage_cnt =
678 k_iter_cnt > max_stage_cnt ? max_stage_cnt : k_iter_cnt;
679 int cta_m_cnt = gemm_m / tile_m;
680 int cta_n_cnt = (gemm_n + tile_n - 1) / tile_n;
681 constexpr int barrier_bytes = (stage_cnt * 16 + 1023) / 1024 * 1024;
682 constexpr int smem_bytes =
683 ((tile_m * 2 + tile_n * 2) * tile_k * stage_cnt + barrier_bytes + 1023) /
684 1024 * 1024;
685
686 dim3 grid(cta_m_cnt, cta_n_cnt, 1);
687 dim3 block_size(256);
688 cudaLaunchConfig_t config;
689 config.gridDim = grid;
690 config.blockDim = block_size;
691 config.dynamicSmemBytes = smem_bytes;
692 config.stream = stream;
693 cudaLaunchAttribute attrs[1];
694 attrs[0].id = cudaLaunchAttributeProgrammaticStreamSerialization;
695 attrs[0].val.programmaticStreamSerializationAllowed =
696 enable_pdl || getEnvEnablePDL();
697 config.numAttrs = 1;
698 config.attrs = attrs;
699 if (smem_bytes >= (48 * 1024)) {
700 cudaFuncSetAttribute(fused_a_gemm_kernel<batch_size, gemm_m, gemm_k, tile_m,
701 tile_n, tile_k, stage_cnt>,
702 cudaFuncAttributeMaxDynamicSharedMemorySize,
703 smem_bytes);
704 }
705 cudaLaunchKernelEx(&config,
706 fused_a_gemm_kernel<batch_size, gemm_m, gemm_k, tile_m,
707 tile_n, tile_k, stage_cnt>,
708 output, mat_a, mat_b, gemm_n, input_row_stride);
709 }
710
711 template <typename T, int kHdIn, int kHdOut, int kTileK = 256, int kTileM = 16>
712 void invokeFusedAGemmForTokens(T* output, T const* mat_a, T const* mat_b,
713 int num_tokens, int input_row_stride,
714 cudaStream_t const stream, bool enable_pdl) {
715 if (num_tokens <= 8) {
716 invokeFusedAGemm<T, kHdIn, kHdOut, 8, kTileK, kTileM>(
717 output, mat_a, mat_b, num_tokens, input_row_stride, stream, enable_pdl);
718 } else {
719 invokeFusedAGemm<T, kHdIn, kHdOut, 16, kTileK, kTileM>(
720 output, mat_a, mat_b, num_tokens, input_row_stride, stream, enable_pdl);
721 }
722 }
723
724 void dsv3_fused_a_gemm(torch::stable::Tensor& output,
725 torch::stable::Tensor const& mat_a,
726 torch::stable::Tensor const& mat_b, bool enable_pdl) {
727 STD_TORCH_CHECK(mat_a.dim() == 2 && mat_b.dim() == 2 && output.dim() == 2);
728 int const num_tokens = mat_a.size(0);
729 int const hd_in = mat_a.size(1);
730 int const hd_out = mat_b.size(1);
731
732 STD_TORCH_CHECK(num_tokens >= 1 && num_tokens <= 16,
733 "required 1 <= mat_a.shape[0] <= 16");
734 STD_TORCH_CHECK(output.size(0) == num_tokens,
735 "required output.shape[0] == mat_a.shape[0]");
736 STD_TORCH_CHECK(output.size(1) == hd_out,
737 "required output.shape[1] == mat_b.shape[1]");
738 STD_TORCH_CHECK(mat_b.size(0) == hd_in,
739 "required mat_b.shape[0] == mat_a.shape[1]");
740
741 STD_TORCH_CHECK(mat_a.get_device_index() == mat_b.get_device_index() &&
742 mat_a.get_device_index() == output.get_device_index(),
743 "mat_a, mat_b, and output must be on the same device");
744
745 STD_TORCH_CHECK(
746 mat_a.stride(1) == 1 && (num_tokens == 1 || mat_a.stride(0) % 8 == 0),
747 "mat_a must have unit inner stride and 16-byte aligned rows");
748 STD_TORCH_CHECK(
749 output.stride(0) == hd_out && output.stride(1) == 1,
750 "output must be a packed row-major [num_tokens, hd_out] tensor");
751 STD_TORCH_CHECK(mat_b.stride(0) == 1 && mat_b.stride(1) == hd_in,
752 "mat_b must be a packed column-major [hd_in, hd_out] tensor");
753
754 STD_TORCH_CHECK(
755 mat_a.scalar_type() == torch::headeronly::ScalarType::BFloat16 &&
756 mat_b.scalar_type() == torch::headeronly::ScalarType::BFloat16,
757 "Only BFloat16 input dtype is supported");
758 STD_TORCH_CHECK(
759 output.scalar_type() == torch::headeronly::ScalarType::BFloat16,
760 "Only BFloat16 output dtype is supported");
761
762 const torch::stable::accelerator::DeviceGuard device_guard(
763 mat_a.get_device_index());
764 STD_TORCH_CHECK(getSMVersion() >= 90, "required CUDA ARCH >= SM_90");
765
766 auto stream = get_current_cuda_stream(mat_a.get_device_index());
767 auto* output_ptr =
768 reinterpret_cast<__nv_bfloat16*>(output.mutable_data_ptr());
769 auto const* mat_a_ptr =
770 reinterpret_cast<__nv_bfloat16 const*>(mat_a.data_ptr());
771 auto const* mat_b_ptr =
772 reinterpret_cast<__nv_bfloat16 const*>(mat_b.data_ptr());
773 int const input_row_stride = mat_a.stride(0);
774
775 #define DISPATCH_DSV3_SHAPE(HD_IN, HD_OUT) \
776 if (hd_in == HD_IN && hd_out == HD_OUT) { \
777 invokeFusedAGemmForTokens<__nv_bfloat16, HD_IN, HD_OUT>( \
778 output_ptr, mat_a_ptr, mat_b_ptr, num_tokens, input_row_stride, \
779 stream, enable_pdl); \
780 return; \
781 }
782
783 // Shapes the Kimi-K3 selector routes to dsv3_fused_a (see the dsv3 winners
784 // in KIMI_K3_PROJECTIONS) plus the DeepSeek V2/V3 QKV A-projection.
785 DISPATCH_DSV3_SHAPE(7168, 1536)
786 DISPATCH_DSV3_SHAPE(7168, 2112)
787 DISPATCH_DSV3_SHAPE(1536, 2304)
788 DISPATCH_DSV3_SHAPE(1536, 4608)
789 DISPATCH_DSV3_SHAPE(7168, 3584)
790 DISPATCH_DSV3_SHAPE(768, 7168)
791 // TP16 dsv3 winners, as (hd_in=K, hd_out=N). TP16 dense down_proj is absent
792 // because hd_in=2112 is not a multiple of any supported tile_k.
793 DISPATCH_DSV3_SHAPE(1536, 1152)
794 DISPATCH_DSV3_SHAPE(7168, 768)
795 DISPATCH_DSV3_SHAPE(7168, 3216)
796 DISPATCH_DSV3_SHAPE(7168, 4224)
797
798 if (hd_in == 6144 && hd_out == 2624) {
799 invokeFusedAGemmForTokens<__nv_bfloat16, 6144, 2624, 256, 32>(
800 output_ptr, mat_a_ptr, mat_b_ptr, num_tokens, input_row_stride, stream,
801 enable_pdl);
802 return;
803 }
804 DISPATCH_DSV3_SHAPE(2048, 2048)
805
806 #ifdef VLLM_K3_BENCH_SHAPES
807 // The selector routes these shapes to CuTe or the default GEMM, so they are
808 // never reached in production. They are compiled only for offline
809 // DSV3-vs-CuTe benchmarking.
810 DISPATCH_DSV3_SHAPE(7168, 6288)
811 DISPATCH_DSV3_SHAPE(1536, 7168)
812 DISPATCH_DSV3_SHAPE(3584, 7168)
813 DISPATCH_DSV3_SHAPE(7168, 8448)
814 DISPATCH_DSV3_SHAPE(7168, 20480)
815 DISPATCH_DSV3_SHAPE(7168, 3072)
816 DISPATCH_DSV3_SHAPE(7168, 12448)
817 DISPATCH_DSV3_SHAPE(3072, 7168)
818 DISPATCH_DSV3_SHAPE(8448, 7168)
819 DISPATCH_DSV3_SHAPE(7168, 16896)
820 DISPATCH_DSV3_SHAPE(7168, 40960)
821 #endif
822
823 #undef DISPATCH_DSV3_SHAPE
824
825 if (hd_in == 128 && hd_out == 1536) {
826 invokeFusedAGemmForTokens<__nv_bfloat16, 128, 1536, 128>(
827 output_ptr, mat_a_ptr, mat_b_ptr, num_tokens, input_row_stride, stream,
828 enable_pdl);
829 return;
830 }
831 if (hd_in == 128 && hd_out == 3072) {
832 invokeFusedAGemmForTokens<__nv_bfloat16, 128, 3072, 128>(
833 output_ptr, mat_a_ptr, mat_b_ptr, num_tokens, input_row_stride, stream,
834 enable_pdl);
835 return;
836 }
837 // TP16 KDA f_b_proj and shared_expert down_proj. Neither hd_in is a multiple
838 // of 256, so both need the 128 tile_k.
839 if (hd_in == 128 && hd_out == 768) {
840 invokeFusedAGemmForTokens<__nv_bfloat16, 128, 768, 128>(
841 output_ptr, mat_a_ptr, mat_b_ptr, num_tokens, input_row_stride, stream,
842 enable_pdl);
843 return;
844 }
845 if (hd_in == 384 && hd_out == 7168) {
846 invokeFusedAGemmForTokens<__nv_bfloat16, 384, 7168, 128>(
847 output_ptr, mat_a_ptr, mat_b_ptr, num_tokens, input_row_stride, stream,
848 enable_pdl);
849 return;
850 }
851 #ifdef VLLM_K3_BENCH_SHAPES
852 if (hd_in == 4224 && hd_out == 7168) {
853 invokeFusedAGemmForTokens<__nv_bfloat16, 4224, 7168, 128>(
854 output_ptr, mat_a_ptr, mat_b_ptr, num_tokens, input_row_stride, stream,
855 enable_pdl);
856 return;
857 }
858 #endif
859
860 STD_TORCH_CHECK(false, "unsupported DSV3 fused-A GEMM shape");
861 }
862
863 STABLE_TORCH_LIBRARY_IMPL(_C, CUDA, m) {
864 m.impl("dsv3_fused_a_gemm", TORCH_BOX(&dsv3_fused_a_gemm));
865 }