main
cpp 546 lines 17.3 KB
Raw
1 // Adapted from
2 // https://github.com/sgl-project/sglang/tree/main/sgl-kernel/csrc/cpu
3
4 // clang-format off
5
6 #include "common.h"
7 #include "gemm.h"
8 #include "vec.h"
9
10 namespace {
11
12 template <typename scalar_t, bool has_bias, int BLOCK_N>
13 struct scale_C {
14 static inline void apply(
15 scalar_t* __restrict__ C,
16 const int32_t* __restrict__ Ctmp,
17 const int32_t* __restrict__ Bcomp,
18 const float* __restrict__ bias,
19 float As,
20 const float* __restrict__ Bs) {
21 TORCH_CHECK(false, "scale_C: scalar path not implemented!");
22 }
23 };
24
25 #if defined(CPU_CAPABILITY_AVX512)
26 template <bool has_bias, int BLOCK_N>
27 struct scale_C<at::BFloat16, has_bias, BLOCK_N> {
28 static inline void apply(
29 at::BFloat16* __restrict__ C,
30 const int32_t* __restrict__ Ctmp,
31 const int32_t* __restrict__ Bcomp,
32 const float* __restrict__ bias,
33 float As,
34 const float* __restrict__ Bs) {
35 constexpr int COLS = BLOCK_N / 16;
36 static_assert(COLS % 2 == 0);
37
38 __m512 vc[COLS];
39 __m512 vd0 = _mm512_set1_ps(As);
40
41 auto compute = [&](auto col) {
42 __m512 vd1 = _mm512_loadu_ps(Bs + col * 16);
43 __m512i vcomp = _mm512_loadu_si512(Bcomp + col * 16);
44 __m512i vc32 = _mm512_loadu_si512(Ctmp + col * 16);
45 vc[col] = _mm512_cvtepi32_ps(_mm512_sub_epi32(vc32, vcomp));
46 if constexpr (has_bias) {
47 __m512 vbias = _mm512_loadu_ps(bias + col * 16);
48 vc[col] = _mm512_fmadd_ps(_mm512_mul_ps(vc[col], vd0), vd1, vbias);
49 } else {
50 vc[col] = _mm512_mul_ps(_mm512_mul_ps(vc[col], vd0), vd1);
51 }
52 };
53 Unroll<COLS>{}(compute);
54
55 auto storec = [&](auto col) {
56 // for COLS = 2, 4 use 512bit store
57 if constexpr (col % 2 == 0) {
58 _mm512_storeu_si512(
59 reinterpret_cast<__m512i*>((C + col * 16)), (__m512i)(_mm512_cvtne2ps_pbh(vc[col + 1], vc[col + 0])));
60 }
61 };
62 Unroll<COLS>{}(storec);
63 }
64 };
65 #endif
66
67 template <typename scalar_t, bool has_bias, int BLOCK_M, int BLOCK_N>
68 struct tinygemm_kernel_nn {
69 static inline void apply(
70 const uint8_t* __restrict__ A,
71 const int8_t* __restrict__ B,
72 scalar_t* __restrict__ C,
73 const float* __restrict__ As,
74 const float* __restrict__ Bs,
75 const int32_t* __restrict__ Bcomp,
76 const float* __restrict__ bias,
77 int64_t K,
78 int64_t lda,
79 int64_t ldb,
80 int64_t ldc) {
81 TORCH_CHECK(false, "tinygemm_kernel_nn: scalar path not implemented!");
82 }
83 };
84
85 #if defined(CPU_CAPABILITY_AVX512)
86 template <bool has_bias, int BLOCK_M, int BLOCK_N>
87 struct tinygemm_kernel_nn<at::BFloat16, has_bias, BLOCK_M, BLOCK_N> {
88 static inline void apply(
89 const uint8_t* __restrict__ A,
90 const int8_t* __restrict__ B,
91 at::BFloat16* __restrict__ C,
92 const float* __restrict__ As,
93 const float* __restrict__ Bs,
94 const int32_t* __restrict__ Bcomp,
95 const float* __restrict__ bias,
96 int64_t K,
97 int64_t lda,
98 int64_t ldb,
99 int64_t ldc) {
100 constexpr int ROWS = BLOCK_M;
101 constexpr int COLS = BLOCK_N / 16;
102 static_assert(COLS % 2 == 0);
103
104 // prefetch distance
105 constexpr int PREFETCH_SIZE_K = 0;
106
107 __m512i va;
108 __m512i vb[COLS];
109 __m512i vc[ROWS * COLS];
110 __m512i vcomp[COLS];
111 __m512 vd0;
112 __m512 vd1[COLS];
113
114 // oops! 4x4 spills but we use 4x2
115 __m512 vbias[COLS];
116
117 // [NOTE]: s8s8 igemm compensation in avx512-vnni
118 //
119 // avx512-vnni has no s8s8, so we need to change s8s8 to u8s8 with compensate:
120 //
121 // a * b = (a + 128) * b - 128 * b
122 // s s u s u s
123 //
124 // 1) 128 * b is pre-computed when packing B to vnni formats
125 // 2) a + 128 is fused when dynamically quantize A
126 //
127 auto loadc = [&](auto i) { vc[i] = _mm512_set1_epi32(0); };
128 Unroll<ROWS * COLS>{}(loadc);
129
130 const int64_t K4 = K >> 2;
131 const int64_t lda4 = lda >> 2;
132 const int64_t ldb4 = ldb; // ldb * 4 >> 2;
133 const int32_t* a_ptr = reinterpret_cast<const int32_t*>(A);
134 const int32_t* b_ptr = reinterpret_cast<const int32_t*>(B);
135
136 auto compute = [&](auto i, int64_t k) {
137 constexpr int row = i / COLS;
138 constexpr int col = i % COLS;
139
140 if constexpr (col == 0) {
141 va = _mm512_set1_epi32(a_ptr[row * lda4 + k]);
142 }
143 if constexpr (row == 0) {
144 vb[col] = _mm512_loadu_si512(b_ptr + k * ldb4 + col * 16);
145 if constexpr (PREFETCH_SIZE_K > 0) {
146 _mm_prefetch(b_ptr + (k + PREFETCH_SIZE_K) * ldb4 + col * 16, _MM_HINT_T0);
147 }
148 }
149 vc[i] = _mm512_dpbusd_epi32(vc[i], va, vb[col]);
150 };
151 for (int64_t k = 0; k < K4; ++k) {
152 Unroll<ROWS * COLS>{}(compute, k);
153 }
154
155 auto storec = [&](auto i) {
156 constexpr int row = i / COLS;
157 constexpr int col = i % COLS;
158
159 // load a scale
160 if constexpr (col == 0) {
161 vd0 = _mm512_set1_ps(As[row]);
162 }
163 // load b scale and vcomp per 2 vectors
164 // also load bias if any
165 if constexpr (row == 0) {
166 if constexpr (col % 2 == 0) {
167 vd1[col + 0] = _mm512_loadu_ps(Bs + col * 16);
168 vd1[col + 1] = _mm512_loadu_ps(Bs + col * 16 + 16);
169 vcomp[col + 0] = _mm512_loadu_si512(Bcomp + col * 16);
170 vcomp[col + 1] = _mm512_loadu_si512(Bcomp + col * 16 + 16);
171 if constexpr (has_bias) {
172 vbias[col + 0] = _mm512_loadu_ps(bias + col * 16);
173 vbias[col + 1] = _mm512_loadu_ps(bias + col * 16 + 16);
174 }
175 }
176 }
177
178 // for COLS = 2, 4 use 512bit store
179 if constexpr (col % 2 == 0) {
180 __m512 vc0 = _mm512_cvtepi32_ps(_mm512_sub_epi32(vc[row * COLS + col + 0], vcomp[col + 0]));
181 __m512 vc1 = _mm512_cvtepi32_ps(_mm512_sub_epi32(vc[row * COLS + col + 1], vcomp[col + 1]));
182 if constexpr (has_bias) {
183 vc0 = _mm512_fmadd_ps(_mm512_mul_ps(vc0, vd0), vd1[col + 0], vbias[col + 0]);
184 vc1 = _mm512_fmadd_ps(_mm512_mul_ps(vc1, vd0), vd1[col + 1], vbias[col + 1]);
185 } else {
186 vc0 = _mm512_mul_ps(_mm512_mul_ps(vc0, vd0), vd1[col + 0]);
187 vc1 = _mm512_mul_ps(_mm512_mul_ps(vc1, vd0), vd1[col + 1]);
188 }
189
190 _mm512_storeu_si512(
191 reinterpret_cast<__m512i*>((C + row * ldc + col * 16)), (__m512i)(_mm512_cvtne2ps_pbh(vc1, vc0)));
192 }
193 };
194 Unroll<ROWS * COLS>{}(storec);
195 }
196 };
197 #endif
198
199 #define LAUNCH_TINYGEMM_KERNEL_NN(MB_SIZE, NB_SIZE) \
200 tinygemm_kernel_nn<scalar_t, has_bias, MB_SIZE, NB_SIZE>::apply( \
201 A + mb_start * lda, \
202 B + nb_start * 4, \
203 C + mb_start * ldc + nb_start, \
204 As + mb_start, \
205 Bs + nb_start, \
206 Bcomp + nb_start, \
207 has_bias ? bias + nb_start : nullptr, \
208 K, \
209 lda, \
210 ldb, \
211 ldc);
212
213 template <typename scalar_t, bool has_bias>
214 void tinygemm_kernel(
215 const uint8_t* __restrict__ A,
216 const int8_t* __restrict__ B,
217 scalar_t* __restrict__ C,
218 int32_t* __restrict__ Ctmp,
219 const float* __restrict__ As,
220 const float* __restrict__ Bs,
221 const float* __restrict__ bias,
222 int64_t M,
223 int64_t N,
224 int64_t K,
225 int64_t lda,
226 int64_t ldb,
227 int64_t ldc,
228 bool brg) {
229 // B compensation
230 const int32_t* Bcomp = reinterpret_cast<const int32_t*>(B + block_size_n() * K);
231
232 if (brg) {
233 constexpr int BLOCK_N = block_size_n();
234 at::native::cpublas::brgemm(M, N, K, lda, ldb, BLOCK_N, /* add_C */ false, A, B, Ctmp);
235
236 // apply compensation and scale
237 for (int64_t m = 0; m < M; ++m) {
238 scale_C<scalar_t, has_bias, BLOCK_N>::apply(C + m * ldc, Ctmp + m * BLOCK_N, Bcomp, bias, As[m], Bs);
239 }
240 return;
241 }
242
243 // pattern: 1-4-16
244 constexpr int64_t BLOCK_M = 4;
245 constexpr int64_t BLOCK_N = 64;
246 const int64_t MB = div_up(M, BLOCK_M);
247 const int64_t NB = div_up(N, BLOCK_N);
248 for (int64_t mb = 0; mb < MB; ++mb) {
249 int64_t mb_start = mb * BLOCK_M;
250 int64_t mb_size = std::min(BLOCK_M, M - mb_start);
251 for (int64_t nb = 0; nb < NB; ++nb) {
252 int64_t nb_start = nb * BLOCK_N;
253 int64_t nb_size = std::min(BLOCK_N, N - nb_start);
254
255 switch (mb_size << 4 | nb_size >> 4) {
256 // mb_size = 1
257 case 0x12:
258 LAUNCH_TINYGEMM_KERNEL_NN(1, 32);
259 break;
260 case 0x14:
261 LAUNCH_TINYGEMM_KERNEL_NN(1, 64);
262 break;
263 // mb_size = 2
264 case 0x22:
265 LAUNCH_TINYGEMM_KERNEL_NN(2, 32);
266 break;
267 case 0x24:
268 LAUNCH_TINYGEMM_KERNEL_NN(2, 64);
269 break;
270 // mb_size = 3
271 case 0x32:
272 LAUNCH_TINYGEMM_KERNEL_NN(3, 32);
273 break;
274 case 0x34:
275 LAUNCH_TINYGEMM_KERNEL_NN(3, 64);
276 break;
277 // mb_size = 4
278 case 0x42:
279 LAUNCH_TINYGEMM_KERNEL_NN(4, 32);
280 break;
281 case 0x44:
282 LAUNCH_TINYGEMM_KERNEL_NN(4, 64);
283 break;
284 default:
285 TORCH_CHECK(false, "Unexpected block size, ", mb_size, "x", "nb_size");
286 }
287 }
288 }
289 }
290
291 template <typename scalar_t>
292 void int8_scaled_mm_kernel_impl(
293 scalar_t* __restrict__ out,
294 const uint8_t* __restrict__ mat1,
295 const int8_t* __restrict__ mat2,
296 const float* __restrict__ scales1,
297 const float* __restrict__ scales2,
298 const float* __restrict__ bias,
299 int64_t M,
300 int64_t N,
301 int64_t K) {
302 constexpr int64_t BLOCK_M = block_size_m();
303 constexpr int64_t BLOCK_N = block_size_n();
304 const int64_t MB = div_up(M, BLOCK_M);
305 const int64_t NB = div_up(N, BLOCK_N);
306
307 const bool use_brgemm = can_use_brgemm<int8_t>(M);
308
309 // K + 4 after compensation
310 const int64_t packed_row_size = get_row_size<int8_t>(K);
311
312 AT_DISPATCH_BOOL(bias != nullptr, has_bias, [&] {
313 parallel_2d(MB, NB, [&](int64_t mb0, int64_t mb1, int64_t nb0, int64_t nb1) {
314 // for brgemm, use int32_t for accumulate
315 alignas(64) int32_t Ctmp[BLOCK_M * BLOCK_N];
316
317 loop_2d<int8_t>(mb0, mb1, nb0, nb1, BLOCK_N * K, [&](int64_t mb, int64_t nb, int64_t nb_offset) {
318 int mb_start = mb * BLOCK_M;
319 int mb_size = std::min(M - mb_start, BLOCK_M);
320 int nb_start = nb * BLOCK_N;
321 int nb_size = std::min(N - nb_start, BLOCK_N);
322
323 tinygemm_kernel<scalar_t, has_bias>(
324 /* A */ mat1 + mb_start * K,
325 /* B */ mat2 + nb_start * packed_row_size /* nb * BLOCK_N * (K + 4) */,
326 /* C */ out + mb_start * N + nb_start,
327 /* Ctmp*/ Ctmp,
328 /* As */ scales1 + mb_start,
329 /* Bs */ scales2 + nb_start,
330 /* bias*/ bias + nb_start,
331 /* M */ mb_size,
332 /* N */ nb_size,
333 /* K */ K,
334 /* lda */ K,
335 /* ldb */ nb_size,
336 /* ldc */ N,
337 /* brg */ use_brgemm);
338 });
339
340 if (use_brgemm) {
341 at::native::cpublas::brgemm_release();
342 }
343 });
344 });
345 }
346
347 } // anonymous namespace
348
349 // tinygemm interface
350 template <typename scalar_t>
351 void tinygemm_kernel(
352 const uint8_t* __restrict__ A,
353 const int8_t* __restrict__ B,
354 scalar_t* __restrict__ C,
355 int32_t* __restrict__ Ctmp,
356 const float* __restrict__ As,
357 const float* __restrict__ Bs,
358 int64_t M,
359 int64_t N,
360 int64_t K,
361 int64_t lda,
362 int64_t ldb,
363 int64_t ldc,
364 bool brg) {
365 tinygemm_kernel<scalar_t, false>(A, B, C, Ctmp, As, Bs, nullptr, M, N, K, lda, ldb, ldc, brg);
366 }
367
368 #define INSTANTIATE_TINYGEMM_TEMPLATE(TYPE) \
369 template void tinygemm_kernel<TYPE>( \
370 const uint8_t* __restrict__ A, \
371 const int8_t* __restrict__ B, \
372 TYPE* __restrict__ C, \
373 int32_t* __restrict__ Ctmp, \
374 const float* __restrict__ As, \
375 const float* __restrict__ Bs, \
376 int64_t M, \
377 int64_t N, \
378 int64_t K, \
379 int64_t lda, \
380 int64_t ldb, \
381 int64_t ldc, \
382 bool brg)
383
384 INSTANTIATE_TINYGEMM_TEMPLATE(at::BFloat16);
385 INSTANTIATE_TINYGEMM_TEMPLATE(at::Half);
386
387 std::tuple<at::Tensor, at::Tensor> per_token_quant_int8_cpu(at::Tensor& A) {
388 CHECK_LAST_DIM_CONTIGUOUS_INPUT(A);
389 CHECK_DIM(2, A);
390
391 int64_t M = A.size(0);
392 int64_t K = A.size(1);
393 int64_t lda = A.stride(0);
394
395 const auto st = A.scalar_type();
396 TORCH_CHECK(st == at::kBFloat16 || st == at::kHalf, "per_token_quant_int8: expect A to be bfloat16 or half.");
397
398 auto Aq = at::empty({M, K}, A.options().dtype(at::kByte));
399 auto As = at::empty({M}, A.options().dtype(at::kFloat));
400
401 AT_DISPATCH_REDUCED_FLOATING_TYPES(st, "per_token_quant_int8", [&] {
402 uint8_t* __restrict__ Aq_data = Aq.data_ptr<uint8_t>();
403 float* __restrict__ As_data = As.data_ptr<float>();
404 const scalar_t* __restrict__ A_data = A.data_ptr<scalar_t>();
405
406 at::parallel_for(0, M, 0, [&](int64_t begin, int64_t end) {
407 for (int64_t m = begin; m < end; ++m) {
408 quantize_row_int8<scalar_t>(Aq_data + m * K, As_data[m], A_data + m * lda, K);
409 }
410 });
411 });
412 return std::make_tuple(Aq, As);
413 }
414
415 // weight : static, per-channel, symmetric
416 // activation : dynamic, per-token, symmetric
417 //
418 // mat1 : [M, K]
419 // mat2 : [N, K]
420 // scales1 : [M]
421 // scales2 : [N]
422 // bias : [N]
423 // out : [M, N]
424 //
425 at::Tensor int8_scaled_mm_cpu(
426 at::Tensor& mat1,
427 at::Tensor& mat2,
428 at::Tensor& scales1,
429 at::Tensor& scales2,
430 const std::optional<at::Tensor>& bias,
431 at::ScalarType out_dtype,
432 bool is_vnni) {
433 auto packed_w = is_vnni ? mat2 : convert_weight_packed(mat2);
434
435 CHECK_INPUT(mat1);
436 CHECK_INPUT(mat2);
437 CHECK_INPUT(scales1);
438 CHECK_INPUT(scales2);
439 CHECK_DIM(2, mat1);
440 CHECK_DIM(2, mat2);
441
442 int64_t M = mat1.size(0);
443 int64_t N = mat2.size(0);
444 int64_t K = mat1.size(1);
445
446 // see [NOTE]: s8s8 igemm compensation in avx512-vnni
447 CHECK_EQ(mat2.size(1), (int64_t)(is_vnni ? K + sizeof(int32_t) : K));
448 CHECK_EQ(scales1.numel(), M);
449 CHECK_EQ(scales2.numel(), N);
450
451 TORCH_CHECK(mat1.scalar_type() == at::kByte, "int8_scaled_mm: expect mat1 to be uint8.");
452 TORCH_CHECK(mat2.scalar_type() == at::kChar, "int8_scaled_mm: expect mat2 to be int8.");
453 TORCH_CHECK(
454 scales1.scalar_type() == at::kFloat && scales2.scalar_type() == at::kFloat,
455 "int8_scaled_mm: expect scales to be float32.");
456
457 auto out = at::empty({M, N}, mat1.options().dtype(out_dtype));
458
459 const bool has_bias = bias.has_value();
460 const float* bias_data = nullptr;
461 if (has_bias) {
462 CHECK_EQ(bias.value().size(0), N);
463 bias_data = bias.value().data_ptr<float>();
464 }
465
466 AT_DISPATCH_REDUCED_FLOATING_TYPES(out_dtype, "int8_scaled_mm_kernel_impl", [&] {
467 int8_scaled_mm_kernel_impl<scalar_t>(
468 out.data_ptr<scalar_t>(),
469 mat1.data_ptr<uint8_t>(),
470 packed_w.data_ptr<int8_t>(),
471 scales1.data_ptr<float>(),
472 scales2.data_ptr<float>(),
473 bias_data,
474 M,
475 N,
476 K);
477 });
478 return out;
479 }
480
481 // fused `per_token_quant_int8_cpu` and `int8_scaled_mm_cpu`
482 at::Tensor int8_scaled_mm_with_quant(
483 at::Tensor& mat1,
484 at::Tensor& mat2,
485 at::Tensor& scales2,
486 const std::optional<at::Tensor>& bias,
487 at::ScalarType out_dtype,
488 bool is_vnni) {
489 auto packed_w = is_vnni ? mat2 : convert_weight_packed(mat2);
490
491 CHECK_LAST_DIM_CONTIGUOUS_INPUT(mat1);
492 CHECK_INPUT(mat2);
493 CHECK_INPUT(scales2);
494 CHECK_DIM(2, mat1);
495 CHECK_DIM(2, mat2);
496
497 int64_t M = mat1.size(0);
498 int64_t N = mat2.size(0);
499 int64_t K = mat1.size(1);
500 int64_t lda = mat1.stride(0);
501
502 // see [NOTE]: s8s8 igemm compensation in avx512-vnni
503 CHECK_EQ(mat2.size(1), (int64_t)(is_vnni ? K + sizeof(int32_t) : K));
504 CHECK_EQ(scales2.numel(), N);
505
506 const auto st = mat1.scalar_type();
507 TORCH_CHECK(st == at::kBFloat16 || st == at::kHalf, "int8_scaled_mm_with_quant: expect A to be bfloat16 or half.");
508 TORCH_CHECK(st == out_dtype, "int8_scaled_mm_with_quant: expect A has same dtype with out_dtype.");
509 TORCH_CHECK(mat2.scalar_type() == at::kChar, "int8_scaled_mm_with_quant: expect mat2 to be int8.");
510 TORCH_CHECK(scales2.scalar_type() == at::kFloat, "int8_scaled_mm_with_quant: expect scales to be float32.");
511
512 const int64_t buffer_size = M * K + M * sizeof(float);
513 auto buffer = at::empty({buffer_size}, mat1.options().dtype(at::kByte));
514 auto out = at::empty({M, N}, mat1.options().dtype(out_dtype));
515
516 const bool has_bias = bias.has_value();
517 const float* bias_data = nullptr;
518 if (has_bias) {
519 CHECK_EQ(bias.value().size(0), N);
520 bias_data = bias.value().data_ptr<float>();
521 }
522
523 AT_DISPATCH_REDUCED_FLOATING_TYPES(out_dtype, "int8_scaled_mm_with_quant_kernel_impl", [&] {
524 uint8_t* __restrict__ Aq_data = buffer.data_ptr<uint8_t>();
525 float* __restrict__ As_data = (float*)((void*)(Aq_data + M * K));
526 const scalar_t* __restrict__ A_data = mat1.data_ptr<scalar_t>();
527
528 at::parallel_for(0, M, 0, [&](int64_t begin, int64_t end) {
529 for (int64_t m = begin; m < end; ++m) {
530 quantize_row_int8<scalar_t>(Aq_data + m * K, As_data[m], A_data + m * lda, K);
531 }
532 });
533
534 int8_scaled_mm_kernel_impl<scalar_t>(
535 out.data_ptr<scalar_t>(),
536 Aq_data,
537 packed_w.data_ptr<int8_t>(),
538 As_data,
539 scales2.data_ptr<float>(),
540 bias_data,
541 M,
542 N,
543 K);
544 });
545 return out;
546 }