| 1 | // Adapted from |
| 2 | // https://github.com/sgl-project/sglang/tree/main/sgl-kernel/csrc/cpu |
| 3 | |
| 4 | // clang-format off |
| 5 | |
| 6 | #include "gemm.h" |
| 7 | |
| 8 | #include "common.h" |
| 9 | #include "vec.h" |
| 10 | |
| 11 | namespace { |
| 12 | |
| 13 | // packed layout: |
| 14 | // quants {N, K} int8_t |
| 15 | // comp {N} int32_t |
| 16 | template <int BLOCK_N> |
| 17 | inline void s8s8_compensation(int8_t* __restrict__ packed, int K) { |
| 18 | #if defined(CPU_CAPABILITY_AVX512) |
| 19 | constexpr int COLS = BLOCK_N / 16; |
| 20 | __m512i vcomp[COLS]; |
| 21 | |
| 22 | for (int col = 0; col < COLS; ++col) { |
| 23 | vcomp[col] = _mm512_setzero_si512(); |
| 24 | } |
| 25 | |
| 26 | const int64_t offset = BLOCK_N * K; |
| 27 | const __m512i off = _mm512_set1_epi8(static_cast<char>(0x80)); |
| 28 | for (int k = 0; k < K / 4; ++k) { |
| 29 | for (int col = 0; col < COLS; ++col) { |
| 30 | __m512i vb = _mm512_loadu_si512((const __m512i*)(packed + k * BLOCK_N * 4 + col * 64)); |
| 31 | vcomp[col] = _mm512_dpbusd_epi32(vcomp[col], off, vb); |
| 32 | } |
| 33 | } |
| 34 | |
| 35 | for (int col = 0; col < COLS; ++col) { |
| 36 | _mm512_storeu_si512((__m512i*)(packed + offset + col * 64), vcomp[col]); |
| 37 | } |
| 38 | #else |
| 39 | TORCH_CHECK(false, "s8s8_compensation not implemented!"); |
| 40 | #endif |
| 41 | } |
| 42 | |
| 43 | // convert to vnni format |
| 44 | // from [N, K] to [K/2, N, 2] for bfloat16 and float16 |
| 45 | template <typename packed_t> |
| 46 | inline void pack_vnni(packed_t* __restrict__ packed, const packed_t* __restrict__ weight, int N, int K) { |
| 47 | const int VNNI_BLK = 2; |
| 48 | for (int n = 0; n < N; ++n) { |
| 49 | for (int k = 0; k < K / VNNI_BLK; ++k) { |
| 50 | for (int d = 0; d < VNNI_BLK; ++d) { |
| 51 | packed[k * N * VNNI_BLK + n * VNNI_BLK + d] = weight[n * K + k * VNNI_BLK + d]; |
| 52 | } |
| 53 | } |
| 54 | } |
| 55 | } |
| 56 | |
| 57 | template <> |
| 58 | inline void pack_vnni<int8_t>(int8_t* __restrict__ packed, const int8_t* __restrict__ weight, int N, int K) { |
| 59 | constexpr int BLOCK_N = block_size_n(); |
| 60 | TORCH_CHECK(N == BLOCK_N); |
| 61 | |
| 62 | const int VNNI_BLK = 4; |
| 63 | for (int n = 0; n < N; ++n) { |
| 64 | for (int k = 0; k < K / VNNI_BLK; ++k) { |
| 65 | for (int d = 0; d < VNNI_BLK; ++d) { |
| 66 | packed[k * N * VNNI_BLK + n * VNNI_BLK + d] = weight[n * K + k * VNNI_BLK + d]; |
| 67 | } |
| 68 | } |
| 69 | } |
| 70 | s8s8_compensation<BLOCK_N>(packed, K); |
| 71 | } |
| 72 | |
| 73 | // uint8_t: mxfp4 or int4 |
| 74 | // pack to vnni2 format as they are computed with bfloat16 |
| 75 | // |
| 76 | // from [N, K'/2, 2] to [K'/2, N, 2], view 2x int4 as unit8: |
| 77 | // from [N, K ] to [K, N ] where K = K'/2 |
| 78 | // |
| 79 | template <> |
| 80 | inline void pack_vnni<uint8_t>(uint8_t* __restrict__ packed, const uint8_t* __restrict__ weight, int N, int K) { |
| 81 | constexpr int BLOCK_N = block_size_n(); |
| 82 | |
| 83 | uint8_t unpacked[2 * BLOCK_N]; |
| 84 | |
| 85 | // 32-way pack (align with BLOCK_N), faster for avx512 unpacking |
| 86 | // |
| 87 | // for a range of (64): |
| 88 | // {0, 1, 2, ..., 63} |
| 89 | // |
| 90 | // original format: |
| 91 | // { 1|0, 3|2, ..., 63|62} |
| 92 | // |
| 93 | // packed format: |
| 94 | // {32|0, 31|1, ..., 63|31} |
| 95 | // |
| 96 | for (int k = 0; k < K; ++k) { |
| 97 | // unpack first |
| 98 | for (int n = 0; n < N; ++n) { |
| 99 | uint8_t value = weight[n * K + k]; |
| 100 | unpacked[n * 2 + 0] = value & 0xF; // lower 4 bits |
| 101 | unpacked[n * 2 + 1] = value >> 4; // higher 4 bits |
| 102 | } |
| 103 | // re-pack to 32-way |
| 104 | for (int n = 0; n < N; ++n) { |
| 105 | packed[k * N + n] = (unpacked[n + BLOCK_N] << 4) | unpacked[n]; |
| 106 | } |
| 107 | } |
| 108 | } |
| 109 | |
| 110 | template <typename scalar_t> |
| 111 | inline void copy_stub(scalar_t* __restrict__ out, const float* __restrict__ input, int64_t size) { |
| 112 | using bVec = at::vec::Vectorized<scalar_t>; |
| 113 | using fVec = at::vec::Vectorized<float>; |
| 114 | constexpr int kVecSize = bVec::size(); |
| 115 | |
| 116 | int64_t d; |
| 117 | #pragma GCC unroll 4 |
| 118 | for (d = 0; d <= size - kVecSize; d += kVecSize) { |
| 119 | auto [data0, data1] = load_float_vec2(input + d); |
| 120 | bVec out_vec = convert_from_float_ext<scalar_t>(data0, data1); |
| 121 | out_vec.store(out + d); |
| 122 | } |
| 123 | for (; d < size; ++d) { |
| 124 | out[d] = static_cast<scalar_t>(input[d]); |
| 125 | } |
| 126 | } |
| 127 | |
| 128 | template <typename scalar_t> |
| 129 | inline void copy_stub(float* __restrict__ out, const scalar_t* __restrict__ input, int64_t size) { |
| 130 | using bVec = at::vec::Vectorized<scalar_t>; |
| 131 | using fVec = at::vec::Vectorized<float>; |
| 132 | constexpr int kVecSize = bVec::size(); |
| 133 | |
| 134 | int64_t d; |
| 135 | #pragma GCC unroll 4 |
| 136 | for (d = 0; d <= size - kVecSize; d += kVecSize) { |
| 137 | auto [data0, data1] = load_float_vec2(input + d); |
| 138 | data0.store(out + d); |
| 139 | data1.store(out + d + fVec::size()); |
| 140 | } |
| 141 | for (; d < size; ++d) { |
| 142 | out[d] = static_cast<float>(input[d]); |
| 143 | } |
| 144 | } |
| 145 | |
| 146 | template <typename scalar_t> |
| 147 | inline void copy_add_stub( |
| 148 | scalar_t* __restrict__ out, const float* __restrict__ input, const float* __restrict__ bias, int64_t size) { |
| 149 | using bVec = at::vec::Vectorized<scalar_t>; |
| 150 | using fVec = at::vec::Vectorized<float>; |
| 151 | constexpr int kVecSize = bVec::size(); |
| 152 | |
| 153 | int64_t d; |
| 154 | #pragma GCC unroll 4 |
| 155 | for (d = 0; d <= size - kVecSize; d += kVecSize) { |
| 156 | auto [data0, data1] = load_float_vec2(input + d); |
| 157 | auto [bias0, bias1] = load_float_vec2(bias + d); |
| 158 | bVec out_vec = convert_from_float_ext<scalar_t>(data0 + bias0, data1 + bias1); |
| 159 | out_vec.store(out + d); |
| 160 | } |
| 161 | for (; d < size; ++d) { |
| 162 | out[d] = static_cast<scalar_t>(input[d] + bias[d]); |
| 163 | } |
| 164 | } |
| 165 | |
| 166 | template <typename scalar_t, bool has_bias> |
| 167 | inline void scalar_sigmoid_and_mul( |
| 168 | scalar_t* __restrict__ out, |
| 169 | const float* __restrict__ input, |
| 170 | const float* __restrict__ bias, |
| 171 | const scalar_t* __restrict__ mul, |
| 172 | int SIZE) { |
| 173 | using bVec = at::vec::Vectorized<scalar_t>; |
| 174 | using fVec = at::vec::Vectorized<float>; |
| 175 | // scalar sigmoid |
| 176 | fVec X; |
| 177 | if constexpr (has_bias) { |
| 178 | assert(bias != nullptr); |
| 179 | X = fVec(input[0] + bias[0]); |
| 180 | } else { |
| 181 | X = fVec(input[0]); |
| 182 | } |
| 183 | X = fast_sigmoid(X); |
| 184 | |
| 185 | // vec mul |
| 186 | constexpr int kVecSize = bVec::size(); |
| 187 | for (int d = 0; d < SIZE; d += kVecSize) { |
| 188 | auto [m_fvec0, m_fvec1] = load_float_vec2(mul + d); |
| 189 | bVec out_vec = convert_from_float_ext<scalar_t>(m_fvec0 * X, m_fvec1 * X); |
| 190 | out_vec.store(out + d); |
| 191 | } |
| 192 | } |
| 193 | |
| 194 | template <typename scalar_t, bool has_bias, int BLOCK_M, int BLOCK_N> |
| 195 | struct tinygemm_kernel_nn { |
| 196 | static inline void apply( |
| 197 | const scalar_t* __restrict__ A, |
| 198 | const scalar_t* __restrict__ B, |
| 199 | scalar_t* __restrict__ C, |
| 200 | const float* __restrict__ bias, |
| 201 | int64_t K, |
| 202 | int64_t lda, |
| 203 | int64_t ldb, |
| 204 | int64_t ldc) { |
| 205 | TORCH_CHECK(false, "tinygemm_kernel_nn: scalar path not implemented!"); |
| 206 | } |
| 207 | }; |
| 208 | |
| 209 | #if defined(CPU_CAPABILITY_AVX512) |
| 210 | template <bool has_bias, int BLOCK_M, int BLOCK_N> |
| 211 | struct tinygemm_kernel_nn<at::BFloat16, has_bias, BLOCK_M, BLOCK_N> { |
| 212 | static inline void apply( |
| 213 | const at::BFloat16* __restrict__ A, |
| 214 | const at::BFloat16* __restrict__ B, |
| 215 | at::BFloat16* __restrict__ C, |
| 216 | const float* __restrict__ bias, |
| 217 | int64_t K, |
| 218 | int64_t lda, |
| 219 | int64_t ldb, |
| 220 | int64_t ldc) { |
| 221 | constexpr int ROWS = BLOCK_M; |
| 222 | constexpr int COLS = BLOCK_N / 16; |
| 223 | |
| 224 | // prefetch distance |
| 225 | constexpr int PREFETCH_SIZE_K = 0; |
| 226 | |
| 227 | __m512bh va; |
| 228 | __m512bh vb[COLS]; |
| 229 | __m512 vc[ROWS * COLS]; |
| 230 | |
| 231 | auto loadc = [&](auto i) { |
| 232 | constexpr int col = i % COLS; |
| 233 | if constexpr (has_bias) { |
| 234 | vc[i] = _mm512_loadu_ps(bias + col * 16); |
| 235 | } else { |
| 236 | vc[i] = _mm512_set1_ps(0.f); |
| 237 | } |
| 238 | }; |
| 239 | Unroll<ROWS * COLS>{}(loadc); |
| 240 | |
| 241 | const int64_t K2 = K >> 1; |
| 242 | const int64_t lda2 = lda >> 1; |
| 243 | const int64_t ldb2 = ldb; // ldb * 2 >> 1; |
| 244 | const float* a_ptr = reinterpret_cast<const float*>(A); |
| 245 | const float* b_ptr = reinterpret_cast<const float*>(B); |
| 246 | |
| 247 | auto compute = [&](auto i, int64_t k) { |
| 248 | constexpr int row = i / COLS; |
| 249 | constexpr int col = i % COLS; |
| 250 | |
| 251 | if constexpr (col == 0) { |
| 252 | va = (__m512bh)(_mm512_set1_ps(a_ptr[row * lda2 + k])); |
| 253 | } |
| 254 | if constexpr (row == 0) { |
| 255 | vb[col] = (__m512bh)(_mm512_loadu_si512(b_ptr + k * ldb2 + col * 16)); |
| 256 | if constexpr (PREFETCH_SIZE_K > 0) { |
| 257 | _mm_prefetch(b_ptr + (k + PREFETCH_SIZE_K) * ldb2 + col * 16, _MM_HINT_T0); |
| 258 | } |
| 259 | } |
| 260 | vc[i] = _mm512_dpbf16_ps(vc[i], va, vb[col]); |
| 261 | }; |
| 262 | for (int64_t k = 0; k < K2; ++k) { |
| 263 | Unroll<ROWS * COLS>{}(compute, k); |
| 264 | } |
| 265 | |
| 266 | auto storec = [&](auto i) { |
| 267 | constexpr int row = i / COLS; |
| 268 | constexpr int col = i % COLS; |
| 269 | // for COLS = 2, 4 use 512bit store |
| 270 | // for COLS = 1, 3 use 256bit store |
| 271 | if constexpr (COLS % 2 == 0) { |
| 272 | if constexpr (col % 2 == 0) { |
| 273 | _mm512_storeu_si512( |
| 274 | reinterpret_cast<__m512i*>((C + row * ldc + col * 16)), |
| 275 | (__m512i)(_mm512_cvtne2ps_pbh(vc[row * COLS + col + 1], vc[row * COLS + col]))); |
| 276 | } |
| 277 | } else { |
| 278 | _mm256_storeu_si256(reinterpret_cast<__m256i*>(C + row * ldc + col * 16), (__m256i)(_mm512_cvtneps_pbh(vc[i]))); |
| 279 | } |
| 280 | }; |
| 281 | Unroll<ROWS * COLS>{}(storec); |
| 282 | } |
| 283 | }; |
| 284 | #endif |
| 285 | |
| 286 | #define LAUNCH_TINYGEMM_KERNEL_NN(MB_SIZE, NB_SIZE) \ |
| 287 | tinygemm_kernel_nn<scalar_t, has_bias, MB_SIZE, NB_SIZE>::apply( \ |
| 288 | A + mb_start * lda, \ |
| 289 | B + nb_start * 2, \ |
| 290 | C + mb_start * ldc + nb_start, \ |
| 291 | has_bias ? bias + nb_start : nullptr, \ |
| 292 | K, \ |
| 293 | lda, \ |
| 294 | ldb, \ |
| 295 | ldc); |
| 296 | |
| 297 | template <typename scalar_t, bool has_bias> |
| 298 | struct brgemm { |
| 299 | static inline void apply( |
| 300 | const scalar_t* __restrict__ A, |
| 301 | const scalar_t* __restrict__ B, |
| 302 | scalar_t* __restrict__ C, |
| 303 | float* __restrict__ Ctmp, |
| 304 | const float* __restrict__ bias, |
| 305 | int64_t M, |
| 306 | int64_t N, |
| 307 | int64_t K, |
| 308 | int64_t lda, |
| 309 | int64_t ldb, |
| 310 | int64_t ldc) { |
| 311 | constexpr int BLOCK_N = block_size_n(); |
| 312 | at::native::cpublas::brgemm(M, N, K, lda, ldb, BLOCK_N, /* add_C */ false, A, B, Ctmp); |
| 313 | |
| 314 | // copy from Ctmp to C |
| 315 | for (int64_t m = 0; m < M; ++m) { |
| 316 | if constexpr (has_bias) { |
| 317 | copy_add_stub(C + m * ldc, Ctmp + m * BLOCK_N, bias, N); |
| 318 | } else { |
| 319 | copy_stub(C + m * ldc, Ctmp + m * BLOCK_N, N); |
| 320 | } |
| 321 | } |
| 322 | } |
| 323 | static inline void apply( |
| 324 | const float* __restrict__ A, |
| 325 | const float* __restrict__ B, |
| 326 | scalar_t* __restrict__ C, |
| 327 | float* __restrict__ Ctmp, |
| 328 | const float* __restrict__ bias, |
| 329 | int64_t M, |
| 330 | int64_t N, |
| 331 | int64_t K, |
| 332 | int64_t lda, |
| 333 | int64_t ldb, |
| 334 | int64_t ldc) { |
| 335 | constexpr int BLOCK_N = block_size_n(); |
| 336 | at::native::cpublas::brgemm(M, N, K, lda, ldb, BLOCK_N, /* add_C */ false, A, B, Ctmp); |
| 337 | } |
| 338 | }; |
| 339 | |
| 340 | template <typename scalar_t, bool has_bias> |
| 341 | void tinygemm_kernel( |
| 342 | const scalar_t* __restrict__ A, |
| 343 | const scalar_t* __restrict__ B, |
| 344 | scalar_t* __restrict__ C, |
| 345 | float* __restrict__ Ctmp, |
| 346 | const float* __restrict__ bias, |
| 347 | int64_t M, |
| 348 | int64_t N, |
| 349 | int64_t K, |
| 350 | int64_t lda, |
| 351 | int64_t ldb, |
| 352 | int64_t ldc, |
| 353 | bool brg) { |
| 354 | if (brg) { |
| 355 | brgemm<scalar_t, has_bias>::apply(A, B, C, Ctmp, bias, M, N, K, lda, ldb, ldc); |
| 356 | return; |
| 357 | } |
| 358 | |
| 359 | // pattern: 1-4-16, N = 16, 32, 48, 64 |
| 360 | constexpr int64_t BLOCK_M = 4; |
| 361 | constexpr int64_t BLOCK_N = 64; |
| 362 | const int64_t MB = div_up(M, BLOCK_M); |
| 363 | const int64_t NB = div_up(N, BLOCK_N); |
| 364 | for (int mb = 0; mb < MB; ++mb) { |
| 365 | int64_t mb_start = mb * BLOCK_M; |
| 366 | int64_t mb_size = std::min(BLOCK_M, M - mb_start); |
| 367 | for (int64_t nb = 0; nb < NB; ++nb) { |
| 368 | int64_t nb_start = nb * BLOCK_N; |
| 369 | int64_t nb_size = std::min(BLOCK_N, N - nb_start); |
| 370 | |
| 371 | switch (mb_size << 4 | nb_size >> 4) { |
| 372 | // mb_size = 1 |
| 373 | case 0x11: |
| 374 | LAUNCH_TINYGEMM_KERNEL_NN(1, 16); |
| 375 | break; |
| 376 | case 0x12: |
| 377 | LAUNCH_TINYGEMM_KERNEL_NN(1, 32); |
| 378 | break; |
| 379 | case 0x13: |
| 380 | LAUNCH_TINYGEMM_KERNEL_NN(1, 48); |
| 381 | break; |
| 382 | case 0x14: |
| 383 | LAUNCH_TINYGEMM_KERNEL_NN(1, 64); |
| 384 | break; |
| 385 | // mb_size = 2 |
| 386 | case 0x21: |
| 387 | LAUNCH_TINYGEMM_KERNEL_NN(2, 16); |
| 388 | break; |
| 389 | case 0x22: |
| 390 | LAUNCH_TINYGEMM_KERNEL_NN(2, 32); |
| 391 | break; |
| 392 | case 0x23: |
| 393 | LAUNCH_TINYGEMM_KERNEL_NN(2, 48); |
| 394 | break; |
| 395 | case 0x24: |
| 396 | LAUNCH_TINYGEMM_KERNEL_NN(2, 64); |
| 397 | break; |
| 398 | // mb_size = 3 |
| 399 | case 0x31: |
| 400 | LAUNCH_TINYGEMM_KERNEL_NN(3, 16); |
| 401 | break; |
| 402 | case 0x32: |
| 403 | LAUNCH_TINYGEMM_KERNEL_NN(3, 32); |
| 404 | break; |
| 405 | case 0x33: |
| 406 | LAUNCH_TINYGEMM_KERNEL_NN(3, 48); |
| 407 | break; |
| 408 | case 0x34: |
| 409 | LAUNCH_TINYGEMM_KERNEL_NN(3, 64); |
| 410 | break; |
| 411 | // mb_size = 4 |
| 412 | case 0x41: |
| 413 | LAUNCH_TINYGEMM_KERNEL_NN(4, 16); |
| 414 | break; |
| 415 | case 0x42: |
| 416 | LAUNCH_TINYGEMM_KERNEL_NN(4, 32); |
| 417 | break; |
| 418 | case 0x43: |
| 419 | LAUNCH_TINYGEMM_KERNEL_NN(4, 48); |
| 420 | break; |
| 421 | case 0x44: |
| 422 | LAUNCH_TINYGEMM_KERNEL_NN(4, 64); |
| 423 | break; |
| 424 | default: |
| 425 | TORCH_CHECK(false, "Unexpected block size, ", mb_size, " x ", nb_size); |
| 426 | } |
| 427 | } |
| 428 | } |
| 429 | } |
| 430 | |
| 431 | template <typename scalar_t, bool has_bias> |
| 432 | void tinygemm_kernel( |
| 433 | const float* __restrict__ A, |
| 434 | const float* __restrict__ B, |
| 435 | scalar_t* __restrict__ C, |
| 436 | float* __restrict__ Ctmp, |
| 437 | const float* __restrict__ bias, |
| 438 | int64_t M, |
| 439 | int64_t N, |
| 440 | int64_t K, |
| 441 | int64_t lda, |
| 442 | int64_t ldb, |
| 443 | int64_t ldc, |
| 444 | bool brg) { |
| 445 | TORCH_CHECK(brg, "Expected to use fp32 brgemm for small N GEMM"); |
| 446 | if (brg) { |
| 447 | brgemm<scalar_t, has_bias>::apply(A, B, C, Ctmp, bias, M, N, K, lda, ldb, ldc); |
| 448 | return; |
| 449 | } |
| 450 | // TODO : add intrinsic path |
| 451 | } |
| 452 | |
| 453 | template <typename scalar_t> |
| 454 | void weight_packed_linear_kernel_impl( |
| 455 | scalar_t* __restrict__ out, |
| 456 | const scalar_t* __restrict__ mat1, |
| 457 | const scalar_t* __restrict__ mat2, |
| 458 | const float* __restrict__ bias, |
| 459 | int64_t M, |
| 460 | int64_t N, |
| 461 | int64_t K, |
| 462 | int64_t mat1_strideM, |
| 463 | int64_t out_strideM) { |
| 464 | constexpr int64_t BLOCK_M = block_size_m(); |
| 465 | constexpr int64_t BLOCK_N = block_size_n(); |
| 466 | const int64_t MB = div_up(M, BLOCK_M); |
| 467 | const int64_t NB = div_up(N, BLOCK_N); |
| 468 | |
| 469 | const bool use_brgemm = can_use_brgemm<scalar_t>(M); |
| 470 | |
| 471 | // parallel on [MB, NB] |
| 472 | AT_DISPATCH_BOOL(bias != nullptr, has_bias, [&] { |
| 473 | parallel_2d(MB, NB, [&](int64_t mb0, int64_t mb1, int64_t nb0, int64_t nb1) { |
| 474 | // for brgemm, use float32 for accumulate |
| 475 | alignas(64) float Ctmp[BLOCK_M * BLOCK_N]; |
| 476 | |
| 477 | loop_2d<scalar_t>(mb0, mb1, nb0, nb1, BLOCK_N * K, [&](int64_t mb, int64_t nb, int64_t nb_offset) { |
| 478 | int64_t mb_start = mb * BLOCK_M; |
| 479 | int64_t mb_size = std::min(M - mb_start, BLOCK_M); |
| 480 | int64_t nb_start = nb * BLOCK_N; |
| 481 | int64_t nb_size = std::min(N - nb_start, BLOCK_N); |
| 482 | |
| 483 | tinygemm_kernel<scalar_t, has_bias>( |
| 484 | /* A */ mat1 + mb_start * mat1_strideM, |
| 485 | /* B */ mat2 + nb_start * K /* nb * BLOCK_N * K */, |
| 486 | /* C */ out + mb_start * out_strideM + nb_start, |
| 487 | /* Ctmp*/ Ctmp, |
| 488 | /* bias*/ bias + nb_start, |
| 489 | /* M */ mb_size, |
| 490 | /* N */ nb_size, |
| 491 | /* K */ K, |
| 492 | /* lda */ mat1_strideM, |
| 493 | /* ldb */ nb_size, |
| 494 | /* ldc */ out_strideM, |
| 495 | /* brg */ use_brgemm); |
| 496 | }); |
| 497 | |
| 498 | if (use_brgemm) { |
| 499 | at::native::cpublas::brgemm_release(); |
| 500 | } |
| 501 | }); |
| 502 | }); |
| 503 | } |
| 504 | |
| 505 | template <typename scalar_t> |
| 506 | void weight_packed_linear_kernel_impl( |
| 507 | scalar_t* __restrict__ out, |
| 508 | const scalar_t* __restrict__ mat1, |
| 509 | const float* __restrict__ mat2, |
| 510 | const float* __restrict__ bias, |
| 511 | const scalar_t* __restrict__ post_mul_mat, |
| 512 | int64_t M, |
| 513 | int64_t N, |
| 514 | int64_t K, |
| 515 | int64_t mat1_strideM, |
| 516 | int64_t out_strideM) { |
| 517 | constexpr int64_t BLOCK_M = block_size_m(); |
| 518 | constexpr int64_t BLOCK_N = block_size_n(); |
| 519 | const int64_t MB = div_up(M, BLOCK_M); |
| 520 | const int64_t NB = div_up(N, BLOCK_N); |
| 521 | |
| 522 | const bool use_brgemm = true; // TODO: add intrinsic path |
| 523 | // parallel on [MB, NB] |
| 524 | AT_DISPATCH_BOOL(bias != nullptr, has_bias, [&] { |
| 525 | parallel_2d(MB, NB, [&](int64_t mb0, int64_t mb1, int64_t nb0, int64_t nb1) { |
| 526 | // for brgemm, use float32 for accumulate |
| 527 | alignas(64) float Atmp[BLOCK_M * K]; |
| 528 | alignas(64) float Ctmp[BLOCK_M * BLOCK_N]; |
| 529 | |
| 530 | loop_2d<float>(mb0, mb1, nb0, nb1, BLOCK_N * K, [&](int64_t mb, int64_t nb, int64_t nb_offset) { |
| 531 | int64_t mb_start = mb * BLOCK_M; |
| 532 | int64_t mb_size = std::min(M - mb_start, BLOCK_M); |
| 533 | int64_t nb_start = nb * BLOCK_N; |
| 534 | int64_t nb_size = std::min(N - nb_start, BLOCK_N); |
| 535 | for (int64_t m = 0; m < mb_size; ++m) { |
| 536 | copy_stub<scalar_t>(Atmp + m * K, mat1 + mb_start * mat1_strideM + m * K, K); |
| 537 | } |
| 538 | tinygemm_kernel<scalar_t, has_bias>( |
| 539 | /* A */ Atmp, |
| 540 | /* B */ mat2 + nb_start * K /* nb * BLOCK_N * K */, |
| 541 | /* C */ out + mb_start * out_strideM + nb_start, |
| 542 | /* Ctmp*/ Ctmp, |
| 543 | /* bias*/ bias + nb_start, |
| 544 | /* M */ mb_size, |
| 545 | /* N */ nb_size, |
| 546 | /* K */ K, |
| 547 | /* lda */ mat1_strideM, |
| 548 | /* ldb */ nb_size, |
| 549 | /* ldc */ out_strideM, |
| 550 | /* brg */ use_brgemm); |
| 551 | |
| 552 | if (post_mul_mat != nullptr) { |
| 553 | for (int64_t m = 0; m < mb_size; ++m) { |
| 554 | scalar_sigmoid_and_mul<scalar_t, has_bias>( |
| 555 | out + mb_start * out_strideM + nb_start + m * out_strideM, |
| 556 | Ctmp + m * BLOCK_N, |
| 557 | bias + nb_start, |
| 558 | post_mul_mat + mb_start * out_strideM + m * out_strideM, |
| 559 | out_strideM); |
| 560 | } |
| 561 | } else { |
| 562 | for (int64_t m = 0; m < mb_size; ++m) { |
| 563 | if constexpr (has_bias) { |
| 564 | copy_add_stub( |
| 565 | out + mb_start * out_strideM + nb_start + m * out_strideM, Ctmp + m * BLOCK_N, bias + nb_start, N); |
| 566 | } else { |
| 567 | copy_stub(out + mb_start * out_strideM + nb_start + m * out_strideM, Ctmp + m * BLOCK_N, N); |
| 568 | } |
| 569 | } |
| 570 | } |
| 571 | }); |
| 572 | |
| 573 | if (use_brgemm) { |
| 574 | at::native::cpublas::brgemm_release(); |
| 575 | } |
| 576 | }); |
| 577 | }); |
| 578 | } |
| 579 | |
| 580 | } // anonymous namespace |
| 581 | |
| 582 | // tinygemm interface |
| 583 | template <typename scalar_t> |
| 584 | void tinygemm_kernel( |
| 585 | const scalar_t* __restrict__ A, |
| 586 | const scalar_t* __restrict__ B, |
| 587 | scalar_t* __restrict__ C, |
| 588 | float* __restrict__ Ctmp, |
| 589 | int64_t M, |
| 590 | int64_t N, |
| 591 | int64_t K, |
| 592 | int64_t lda, |
| 593 | int64_t ldb, |
| 594 | int64_t ldc, |
| 595 | bool brg) { |
| 596 | tinygemm_kernel<scalar_t, false>(A, B, C, Ctmp, nullptr, M, N, K, lda, ldb, ldc, brg); |
| 597 | } |
| 598 | |
| 599 | #define INSTANTIATE_TINYGEMM_TEMPLATE(TYPE) \ |
| 600 | template void tinygemm_kernel<TYPE>( \ |
| 601 | const TYPE* __restrict__ A, \ |
| 602 | const TYPE* __restrict__ B, \ |
| 603 | TYPE* __restrict__ C, \ |
| 604 | float* __restrict__ Ctmp, \ |
| 605 | int64_t M, \ |
| 606 | int64_t N, \ |
| 607 | int64_t K, \ |
| 608 | int64_t lda, \ |
| 609 | int64_t ldb, \ |
| 610 | int64_t ldc, \ |
| 611 | bool brg) |
| 612 | |
| 613 | INSTANTIATE_TINYGEMM_TEMPLATE(at::BFloat16); |
| 614 | INSTANTIATE_TINYGEMM_TEMPLATE(at::Half); |
| 615 | |
| 616 | at::Tensor convert_weight_packed(at::Tensor& weight) { |
| 617 | // for 3d moe weights |
| 618 | // weight : [E, OC, IC] |
| 619 | // w1 : [E, 2N, K] |
| 620 | // w2 : [E, K, N] |
| 621 | CHECK_INPUT(weight); |
| 622 | |
| 623 | const int64_t ndim = weight.ndimension(); |
| 624 | TORCH_CHECK(ndim == 2 || ndim == 3, "expect weight to be 2d or 3d, got ", ndim, "d tensor."); |
| 625 | |
| 626 | if (ndim == 2 && weight.size(0) < TILE_N) { |
| 627 | // for 2D weight and small OC shape, we use fma linear path, which needs transpose not pack |
| 628 | return weight.to(at::kFloat).t().contiguous(); |
| 629 | } |
| 630 | |
| 631 | const auto st = weight.scalar_type(); |
| 632 | const int64_t E = ndim == 3 ? weight.size(0) : 1; |
| 633 | const int64_t OC = ndim == 3 ? weight.size(1) : weight.size(0); |
| 634 | const int64_t IC = ndim == 3 ? weight.size(2) : weight.size(1); |
| 635 | |
| 636 | // mxfp4 or int4 are packed with uint8 |
| 637 | const int64_t actual_IC = st == at::kByte ? IC * 2 : IC; |
| 638 | |
| 639 | // we handle 2 TILE_N at a time. |
| 640 | TORCH_CHECK(OC % TILE_N == 0, "invalid weight out features ", OC); |
| 641 | TORCH_CHECK(actual_IC % TILE_K == 0, "invalid weight input features ", actual_IC); |
| 642 | |
| 643 | constexpr int64_t BLOCK_N = block_size_n(); |
| 644 | const int64_t NB = div_up(OC, BLOCK_N); |
| 645 | |
| 646 | // use phony sizes here [E, OC, IC], for each [E], [OC, IC] -> [IC / 2, OC, 2] |
| 647 | auto packed_weight = at::empty({}, weight.options()); |
| 648 | const int64_t stride = OC * IC; |
| 649 | |
| 650 | // Note: for `kByte` (uint8), it represents either `mxfp4` or `int4`. |
| 651 | TORCH_CHECK( |
| 652 | st == at::kBFloat16 || st == at::kHalf || st == at::kChar || st == at::kFloat8_e4m3fn || st == at::kByte, |
| 653 | "expect weight to be bfloat16, float16, int8, fp8_e4m3 or uint8(mxfp4 or int4)."); |
| 654 | |
| 655 | CPU_DISPATCH_PACKED_TYPES(st, [&] { |
| 656 | // adjust most inner dimension size |
| 657 | const int packed_row_size = get_row_size<packed_t>(actual_IC); |
| 658 | auto sizes = weight.sizes().vec(); |
| 659 | sizes[ndim - 1] = packed_row_size; |
| 660 | packed_weight.resize_(sizes); |
| 661 | |
| 662 | const packed_t* w_data = weight.data_ptr<packed_t>(); |
| 663 | packed_t* packed_data = packed_weight.data_ptr<packed_t>(); |
| 664 | |
| 665 | // parallel on {E, NB} |
| 666 | at::parallel_for(0, E * NB, 0, [&](int64_t begin, int64_t end) { |
| 667 | int64_t e{0}, nb{0}; |
| 668 | data_index_init(begin, e, E, nb, NB); |
| 669 | |
| 670 | for (int64_t i = begin; i < end; ++i) { |
| 671 | UNUSED(i); |
| 672 | |
| 673 | int64_t n = nb * BLOCK_N; |
| 674 | int64_t n_size = std::min(BLOCK_N, OC - n); |
| 675 | pack_vnni<packed_t>( |
| 676 | packed_data + e * OC * packed_row_size + n * packed_row_size, w_data + e * stride + n * IC, n_size, IC); |
| 677 | |
| 678 | // move to the next index |
| 679 | data_index_step(e, E, nb, NB); |
| 680 | } |
| 681 | }); |
| 682 | }); |
| 683 | return packed_weight; |
| 684 | } |
| 685 | |
| 686 | at::Tensor convert_scale_packed(at::Tensor& scale) { |
| 687 | CHECK_INPUT(scale); |
| 688 | |
| 689 | const int64_t ndim = scale.ndimension(); |
| 690 | TORCH_CHECK(ndim == 2 || ndim == 3, "expect scale to be 2d or 3d, got ", ndim, "d tensor."); |
| 691 | const auto st = scale.scalar_type(); |
| 692 | const int64_t E = ndim == 3 ? scale.size(0) : 1; |
| 693 | const int64_t N = ndim == 3 ? scale.size(1) : scale.size(0); |
| 694 | // number of groups, e.g. K/32 |
| 695 | const int64_t G = ndim == 3 ? scale.size(2) : scale.size(1); |
| 696 | |
| 697 | constexpr int64_t BLOCK_N = block_size_n(); |
| 698 | TORCH_CHECK(N % BLOCK_N == 0, "invalid weight out features ", N); |
| 699 | const int64_t NB = N / BLOCK_N; |
| 700 | |
| 701 | auto packed_scale = at::empty_like(scale); |
| 702 | TORCH_CHECK(st == at::kByte, "expect scale to be uint8."); |
| 703 | |
| 704 | const uint8_t* s_data = scale.data_ptr<uint8_t>(); |
| 705 | uint8_t* packed_data = packed_scale.data_ptr<uint8_t>(); |
| 706 | |
| 707 | // parallel on src {E, NB, BLOCK_N, G}, dst {E, NB, G, BLOCK_N} |
| 708 | at::parallel_for(0, E * NB * BLOCK_N * G, 0, [&](int64_t begin, int64_t end) { |
| 709 | int64_t e{0}, nb{0}, n{0}, g{0}; |
| 710 | data_index_init(begin, e, E, nb, NB, n, BLOCK_N, g, G); |
| 711 | |
| 712 | for (int64_t i = begin; i < end; ++i) { |
| 713 | packed_data[e * N * G + nb * G * BLOCK_N + g * BLOCK_N + n] = s_data[i]; |
| 714 | // move to the next index |
| 715 | data_index_step(e, E, nb, NB, n, BLOCK_N, g, G); |
| 716 | } |
| 717 | }); |
| 718 | return packed_scale; |
| 719 | } |
| 720 | |
| 721 | // mat1 : [*, K] |
| 722 | // mat2 : [N, K] ([K, N] if use_fma_gemm) |
| 723 | // bias : [N] |
| 724 | // out : [*, N] |
| 725 | // |
| 726 | void |
| 727 | weight_packed_linear(at::Tensor& out, const at::Tensor& mat1, const at::Tensor& mat2, const std::optional<at::Tensor>& bias, bool is_vnni) { |
| 728 | TORCH_CHECK(is_vnni); |
| 729 | auto packed_w = mat2; |
| 730 | bool use_fma_gemm = false; |
| 731 | if (packed_w.scalar_type() == at::kFloat) { |
| 732 | use_fma_gemm = true; |
| 733 | } |
| 734 | |
| 735 | CHECK_LAST_DIM_CONTIGUOUS_INPUT(mat1); |
| 736 | CHECK_INPUT(mat2); |
| 737 | const int64_t ndim = mat1.ndimension(); |
| 738 | auto input_sizes = mat1.sizes().vec(); |
| 739 | int64_t N = use_fma_gemm ? mat2.size(1) : mat2.size(0); |
| 740 | int64_t K = use_fma_gemm ? mat1.size(1) : mat2.size(1); |
| 741 | int64_t M = use_fma_gemm ? mat1.size(0) : mat1.numel() / K; |
| 742 | CHECK_DIM(2, mat2); |
| 743 | if (use_fma_gemm) { |
| 744 | CHECK_DIM(2, mat1); |
| 745 | } else { |
| 746 | CHECK_EQ(mat1.size(ndim - 1), K); |
| 747 | } |
| 748 | |
| 749 | auto dispatch_type = mat1.scalar_type(); |
| 750 | // strides |
| 751 | int64_t out_strideM = out.stride(0); |
| 752 | int64_t mat1_strideM = mat1.stride(-2); |
| 753 | |
| 754 | const bool has_bias = bias.has_value(); |
| 755 | const float* bias_data = nullptr; |
| 756 | if (has_bias) { |
| 757 | CHECK_EQ(bias.value().size(0), N); |
| 758 | bias_data = bias.value().data_ptr<float>(); |
| 759 | } |
| 760 | |
| 761 | AT_DISPATCH_REDUCED_FLOATING_TYPES(dispatch_type, "weight_packed_linear_kernel_impl", [&] { |
| 762 | if (use_fma_gemm) { |
| 763 | weight_packed_linear_kernel_impl<scalar_t>( |
| 764 | out.data_ptr<scalar_t>(), |
| 765 | mat1.data_ptr<scalar_t>(), |
| 766 | packed_w.data_ptr<float>(), |
| 767 | bias_data, |
| 768 | nullptr, |
| 769 | M, |
| 770 | N, |
| 771 | K, |
| 772 | mat1_strideM, |
| 773 | out_strideM); |
| 774 | } else { |
| 775 | weight_packed_linear_kernel_impl<scalar_t>( |
| 776 | out.data_ptr<scalar_t>(), |
| 777 | mat1.data_ptr<scalar_t>(), |
| 778 | packed_w.data_ptr<scalar_t>(), |
| 779 | bias_data, |
| 780 | M, |
| 781 | N, |
| 782 | K, |
| 783 | mat1_strideM, |
| 784 | out_strideM); |
| 785 | } |
| 786 | }); |
| 787 | } |
| 788 | |
| 789 | // mat1 : [M, K] |
| 790 | // mat2 : [K, 1] |
| 791 | // post_mul_mat : [M, K] |
| 792 | // bias : [N] |
| 793 | // out : [M, N] |
| 794 | // |
| 795 | at::Tensor fused_linear_sigmoid_mul( |
| 796 | at::Tensor& mat1, |
| 797 | at::Tensor& mat2, |
| 798 | const std::optional<at::Tensor>& bias, |
| 799 | bool is_vnni, |
| 800 | const at::Tensor& post_mul_mat) { |
| 801 | auto packed_w = is_vnni ? mat2 : convert_weight_packed(mat2); |
| 802 | TORCH_CHECK(packed_w.scalar_type() == at::kFloat, "fused_linear_sigmoid_mul requires packed float weight") |
| 803 | |
| 804 | int64_t M = mat1.size(0); |
| 805 | int64_t K = mat1.size(1); |
| 806 | int64_t N = mat2.size(1); |
| 807 | |
| 808 | CHECK_LAST_DIM_CONTIGUOUS_INPUT(mat1); |
| 809 | CHECK_INPUT(mat2); |
| 810 | CHECK_DIM(2, mat1); |
| 811 | CHECK_DIM(2, mat2); |
| 812 | |
| 813 | int64_t out_strideM = post_mul_mat.size(1); |
| 814 | int64_t mat1_strideM = mat1.stride(0); |
| 815 | auto dispatch_type = mat1.scalar_type(); |
| 816 | auto out = at::empty({M, out_strideM}, mat1.options()); |
| 817 | |
| 818 | TORCH_CHECK( |
| 819 | N == 1 && out_strideM % 32 == 0, |
| 820 | "post_mul_mat tensor size(1) should be 32 dividable, and the mat2 OC=1 (Mx1 as linear output shape)") |
| 821 | |
| 822 | const bool has_bias = bias.has_value(); |
| 823 | const float* bias_data = nullptr; |
| 824 | if (has_bias) { |
| 825 | CHECK_EQ(bias.value().size(0), N); |
| 826 | bias_data = bias.value().data_ptr<float>(); |
| 827 | } |
| 828 | |
| 829 | AT_DISPATCH_REDUCED_FLOATING_TYPES(dispatch_type, "fused_linear_sigmoid_mul", [&] { |
| 830 | weight_packed_linear_kernel_impl<scalar_t>( |
| 831 | out.data_ptr<scalar_t>(), |
| 832 | mat1.data_ptr<scalar_t>(), |
| 833 | packed_w.data_ptr<float>(), |
| 834 | bias_data, |
| 835 | post_mul_mat.data_ptr<scalar_t>(), |
| 836 | M, |
| 837 | N, |
| 838 | K, |
| 839 | mat1_strideM, |
| 840 | out_strideM); |
| 841 | }); |
| 842 | |
| 843 | return out; |
| 844 | } |