← Back to C-Kernel-Engine Docs Doxygen Source Documentation
 
Loading...
Searching...
No Matches
gemm_kernels_q6k_q8k.c
Go to the documentation of this file.
1/**
2 * @file gemm_kernels_q6k_q8k.c
3 * @brief Q6_K (weights) x Q8_K (activations) kernels for inference
4 *
5 * CK-ENGINE KERNEL RULES:
6 * =======================
7 * 1. NO malloc/free - memory via bump allocator, pointers passed in
8 * 2. NO OpenMP - parallelization at orchestrator/codegen layer
9 * 3. API must define: inputs, outputs, workspace, and memory layouts
10 * 4. Pure computation - deterministic, no side effects
11 *
12 * After changes: make test && make llamacpp-parity-full
13 *
14 * Implements decode-style matvec/matmul where weights are Q6_K and the
15 * activations are quantized on-the-fly to Q8_K. This is inference-only;
16 * no backward pass is provided here.
17 *
18 * Q6_K Format (256 weights per block):
19 * - d: FP16 super-block scale
20 * - ql: 128 bytes (low 4 bits of each weight)
21 * - qh: 64 bytes (high 2 bits of each weight)
22 * - scales: 16 int8 sub-block scales
23 *
24 * Q8_K Format (256 weights per block):
25 * - d: FP32 scale
26 * - qs: 256 int8 values
27 * - bsums: 16 int16 block sums
28 */
29
30#include <assert.h>
31#include <math.h>
32#include <string.h>
33#include <stdint.h>
34#include <stddef.h>
35#include <stdlib.h>
36
37#include "ckernel_engine.h"
38#include "ckernel_quant.h"
39
40/* Include SIMD headers based on available extensions */
41#if defined(__AVX512F__) || defined(__AVX2__) || defined(__AVX__) || defined(__SSE4_1__) || defined(__SSSE3__)
42#include <immintrin.h>
43#endif
44#if defined(__ARM_NEON) || defined(__aarch64__)
45#include <arm_neon.h>
46#endif
47
48typedef struct {
49 ggml_half d;
50 int8_t scales[16];
51 uint8_t qs[QK_K];
52} block_q6_K_prepared;
53
54_Static_assert(sizeof(block_q6_K_prepared) == 274,
55 "Q6_K prepared-size contract changed");
56
58{
59 return sizeof(block_q6_K_prepared);
60}
61
63{
64#if defined(__AVX512F__) && defined(__AVX512BW__) && \
65 defined(__AVX512VNNI__)
66 return "q6_k_prepared_avx512_vnni_exact";
67#elif defined(__AVX2__)
68 return "q6_k_prepared_avx2_exact";
69#else
70 return "q6_k_prepared_unavailable";
71#endif
72}
73
74void ck_q6_k_prepare_weight(const void *src, void *dst, int N, int K)
75{
76 if (!src || !dst || N <= 0 || K <= 0 || (K % QK_K) != 0) return;
77 const block_q6_K *input = (const block_q6_K *)src;
78 block_q6_K_prepared *output = (block_q6_K_prepared *)dst;
79 const size_t blocks = (size_t)N * (size_t)(K / QK_K);
80
81 for (size_t b = 0; b < blocks; ++b) {
82 output[b].d = input[b].d;
83 memcpy(output[b].scales, input[b].scales, sizeof(output[b].scales));
84 for (int n = 0; n < QK_K; n += 128) {
85 const uint8_t *ql = input[b].ql + n / 2;
86 const uint8_t *qh = input[b].qh + n / 4;
87 for (int l = 0; l < 32; ++l) {
88 output[b].qs[n + l + 0] =
89 (uint8_t)((ql[l] & 0x0f) | (((qh[l] >> 0) & 3) << 4));
90 output[b].qs[n + l + 32] =
91 (uint8_t)((ql[l + 32] & 0x0f) | (((qh[l] >> 2) & 3) << 4));
92 output[b].qs[n + l + 64] =
93 (uint8_t)((ql[l] >> 4) | (((qh[l] >> 4) & 3) << 4));
94 output[b].qs[n + l + 96] =
95 (uint8_t)((ql[l + 32] >> 4) | (((qh[l] >> 6) & 3) << 4));
96 }
97 }
98 }
99}
100
101/* Forward declarations for SIMD implementations */
102void gemv_q6_k_q8_k_avx512(float *y, const void *W, const void *x_q8, int M, int K);
103void gemv_q6_k_q8_k_avx512_vbmi(float *y, const void *W, const void *x_q8, int M, int K);
104void gemv_q6_k_q8_k_avx2(float *y, const void *W, const void *x_q8, int M, int K);
105void gemv_q6_k_q8_k_avx(float *y, const void *W, const void *x_q8, int M, int K);
106void gemv_q6_k_q8_k_sse(float *y, const void *W, const void *x_q8, int M, int K);
107
108static int ck_q6k_q8k_force_ref(void)
109{
110 static int cached = -1;
111 if (cached < 0) {
112 const char *env = getenv("CK_DEBUG_Q6K_Q8K_REF");
113 cached = (env && env[0] && env[0] != '0') ? 1 : 0;
114 }
115 return cached;
116}
117
118/* ============================================================================
119 * Reference Implementation
120 * ============================================================================ */
121
122/**
123 * @brief Scalar dot product for Q6_K x Q8_K
124 *
125 * Q6_K layout: 256 weights per block
126 * - ql[0..127]: low 4 bits for all 256 weights (packed 2 per byte)
127 * - qh[0..63]: high 2 bits for all 256 weights (packed 4 per byte)
128 * - scales[0..15]: int8 scale for each 16-weight sub-block
129 * - d: FP16 super-block scale
130 *
131 * The dequantization formula for each weight is:
132 * weight = d * scale[sub] * (q6_value - 32)
133 * where q6_value is the 6-bit unsigned value (0..63).
134 */
135static float dot_q6_k_q8_k_ref(const block_q6_K *w,
136 const block_q8_K *x,
137 int K)
138{
139 const int nb = K / QK_K;
140 float sumf = 0.0f;
141 float sums[8] = {0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
142
143 for (int i = 0; i < nb; ++i) {
144 const float d = GGML_FP16_TO_FP32(w[i].d) * x[i].d;
145
146 const uint8_t *ql = w[i].ql;
147 const uint8_t *qh = w[i].qh;
148 const int8_t *sc = w[i].scales;
149 const int8_t *q8 = x[i].qs;
150 int32_t aux32[8] = {0, 0, 0, 0, 0, 0, 0, 0};
151
152 /* Process 256 weights in 2 iterations of 128 */
153 for (int n = 0; n < QK_K; n += 128) {
154 /* Each iteration processes 128 weights:
155 * - ql[0..63] contains low 4 bits
156 * - qh[0..31] contains high 2 bits
157 * - Interleaved pattern: weights 0-31, 32-63, 64-95, 96-127
158 */
159 for (int l = 0; l < 32; ++l) {
160 /* Sub-block index: each scale covers 16 weights */
161 const int is = l / 16;
162
163 /* Extract 6-bit values from packed format */
164 /* q1: weights l+0 (low nibble of ql[l], bits 0-1 of qh[l]) */
165 const int8_t q1 = (int8_t)((ql[l + 0] & 0xF) | (((qh[l] >> 0) & 3) << 4)) - 32;
166 /* q2: weights l+32 (low nibble of ql[l+32], bits 2-3 of qh[l]) */
167 const int8_t q2 = (int8_t)((ql[l + 32] & 0xF) | (((qh[l] >> 2) & 3) << 4)) - 32;
168 /* q3: weights l+64 (high nibble of ql[l], bits 4-5 of qh[l]) */
169 const int8_t q3 = (int8_t)((ql[l + 0] >> 4) | (((qh[l] >> 4) & 3) << 4)) - 32;
170 /* q4: weights l+96 (high nibble of ql[l+32], bits 6-7 of qh[l]) */
171 const int8_t q4 = (int8_t)((ql[l + 32] >> 4) | (((qh[l] >> 6) & 3) << 4)) - 32;
172
173 aux32[l & 7] += (int)sc[is + 0] * (int)q1 * (int)q8[l + 0];
174 aux32[l & 7] += (int)sc[is + 2] * (int)q2 * (int)q8[l + 32];
175 aux32[l & 7] += (int)sc[is + 4] * (int)q3 * (int)q8[l + 64];
176 aux32[l & 7] += (int)sc[is + 6] * (int)q4 * (int)q8[l + 96];
177 }
178 q8 += 128;
179 ql += 64;
180 qh += 32;
181 sc += 8;
182 }
183
184 for (int l = 0; l < 8; ++l) {
185 sums[l] += d * (float)aux32[l];
186 }
187 }
188
189 for (int l = 0; l < 8; ++l) {
190 sumf += sums[l];
191 }
192 return sumf;
193}
194
195#if defined(__ARM_NEON) || defined(__aarch64__)
196static float dot_q6_k_q8_k_neon(const block_q6_K *w,
197 const block_q8_K *x,
198 int K)
199{
200 const int nb = K / QK_K;
201 float sumf = 0.0f;
202
203 for (int i = 0; i < nb; ++i) {
204 const float d = GGML_FP16_TO_FP32(w[i].d) * x[i].d;
205
206 const uint8_t *ql = w[i].ql;
207 const uint8_t *qh = w[i].qh;
208 const int8_t *sc = w[i].scales;
209 const int8_t *q8 = x[i].qs;
210
211 int8_t wvals[QK_K];
212 int8_t svals[QK_K];
213
214 for (int n = 0; n < QK_K; n += 128) {
215 for (int l = 0; l < 32; ++l) {
216 const int is = l / 16;
217
218 const int8_t q1 = (int8_t)((ql[l + 0] & 0xF) | (((qh[l] >> 0) & 3) << 4)) - 32;
219 const int8_t q2 = (int8_t)((ql[l + 32] & 0xF) | (((qh[l] >> 2) & 3) << 4)) - 32;
220 const int8_t q3 = (int8_t)((ql[l + 0] >> 4) | (((qh[l] >> 4) & 3) << 4)) - 32;
221 const int8_t q4 = (int8_t)((ql[l + 32] >> 4) | (((qh[l] >> 6) & 3) << 4)) - 32;
222
223 const int base = n;
224 wvals[base + l + 0] = q1;
225 wvals[base + l + 32] = q2;
226 wvals[base + l + 64] = q3;
227 wvals[base + l + 96] = q4;
228
229 svals[base + l + 0] = sc[is + 0];
230 svals[base + l + 32] = sc[is + 2];
231 svals[base + l + 64] = sc[is + 4];
232 svals[base + l + 96] = sc[is + 6];
233 }
234
235 ql += 64;
236 qh += 32;
237 sc += 8;
238 }
239
240 int32x4_t acc = vdupq_n_s32(0);
241 for (int j = 0; j < QK_K; j += 16) {
242 const int8x16_t wv = vld1q_s8(&wvals[j]);
243 const int8x16_t sv = vld1q_s8(&svals[j]);
244 const int8x16_t xv = vld1q_s8(&q8[j]);
245
246 const int16x8_t ws0 = vmull_s8(vget_low_s8(wv), vget_low_s8(sv));
247 const int16x8_t ws1 = vmull_s8(vget_high_s8(wv), vget_high_s8(sv));
248 const int16x8_t x0 = vmovl_s8(vget_low_s8(xv));
249 const int16x8_t x1 = vmovl_s8(vget_high_s8(xv));
250
251 int32x4_t p0 = vmull_s16(vget_low_s16(ws0), vget_low_s16(x0));
252 p0 = vmlal_s16(p0, vget_high_s16(ws0), vget_high_s16(x0));
253
254 int32x4_t p1 = vmull_s16(vget_low_s16(ws1), vget_low_s16(x1));
255 p1 = vmlal_s16(p1, vget_high_s16(ws1), vget_high_s16(x1));
256
257 acc = vaddq_s32(acc, p0);
258 acc = vaddq_s32(acc, p1);
259 }
260
261 int32_t lanes[4];
262 vst1q_s32(lanes, acc);
263 sumf += d * (float)(lanes[0] + lanes[1] + lanes[2] + lanes[3]);
264 }
265
266 return sumf;
267}
268#endif
269
270void gemv_q6_k_q8_k_ref(float *y,
271 const void *W,
272 const void *x_q8,
273 int M, int K)
274{
275 if (!y || !W || !x_q8 || M <= 0 || K <= 0) {
276 return;
277 }
278
279 const block_q6_K *blocks = (const block_q6_K *)W;
280 const block_q8_K *x = (const block_q8_K *)x_q8;
281 const int blocks_per_row = K / QK_K;
282
283 for (int row = 0; row < M; ++row) {
284 const block_q6_K *w_row = blocks + (size_t)row * (size_t)blocks_per_row;
285 y[row] = dot_q6_k_q8_k_ref(w_row, x, K);
286 }
287}
288
289/* ============================================================================
290 * SSE4.1 Implementation (for Ivy Bridge and older AVX-without-AVX2 CPUs)
291 *
292 * Uses 128-bit SSE operations with maddubs for integer multiply-add.
293 * Handles the -32 offset using bsums from Q8_K.
294 * ============================================================================ */
295
296#if defined(__SSSE3__)
297
298/* Scale shuffle indices for Q6_K - maps scale index to 16-byte shuffle pattern */
299static const int8_t q6k_scale_shuffle[8][16] = {
300 { 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1 }, /* is=0: scales[0,1] */
301 { 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3 }, /* is=1: scales[2,3] */
302 { 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5 }, /* is=2: scales[4,5] */
303 { 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7 }, /* is=3: scales[6,7] */
304 { 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9 }, /* is=4: scales[8,9] */
305 {10,10,10,10,10,10,10,10,11,11,11,11,11,11,11,11 }, /* is=5: scales[10,11] */
306 {12,12,12,12,12,12,12,12,13,13,13,13,13,13,13,13 }, /* is=6: scales[12,13] */
307 {14,14,14,14,14,14,14,14,15,15,15,15,15,15,15,15 }, /* is=7: scales[14,15] */
308};
309
310static float dot_q6_k_q8_k_sse(const block_q6_K *w,
311 const block_q8_K *x,
312 int K)
313{
314 const int nb = K / QK_K;
315 const __m128i m3 = _mm_set1_epi8(3);
316 const __m128i m15 = _mm_set1_epi8(15);
317
318 __m128 acc = _mm_setzero_ps();
319
320 for (int i = 0; i < nb; ++i) {
321 const float d = GGML_FP16_TO_FP32(w[i].d) * x[i].d;
322
323 const uint8_t *ql = w[i].ql;
324 const uint8_t *qh = w[i].qh;
325 const int8_t *q8 = x[i].qs;
326
327 /* Load scales and precompute the -32 offset contribution using bsums */
328 const __m128i scales = _mm_loadu_si128((const __m128i *)w[i].scales);
329 const __m128i q8sums_0 = _mm_loadu_si128((const __m128i *)x[i].bsums);
330 const __m128i q8sums_1 = _mm_loadu_si128((const __m128i *)x[i].bsums + 1);
331
332 /* Compute: sum(scale * bsum) * 32 for the -32 offset */
333 const __m128i scales_16_0 = _mm_cvtepi8_epi16(scales);
334 const __m128i scales_16_1 = _mm_cvtepi8_epi16(_mm_bsrli_si128(scales, 8));
335 const __m128i q8sclsub_0 = _mm_slli_epi32(_mm_madd_epi16(q8sums_0, scales_16_0), 5);
336 const __m128i q8sclsub_1 = _mm_slli_epi32(_mm_madd_epi16(q8sums_1, scales_16_1), 5);
337
338 __m128i sumi_0 = _mm_setzero_si128();
339 __m128i sumi_1 = _mm_setzero_si128();
340
341 int is = 0;
342
343 /* Process 256 weights in 2 iterations of 128 */
344 for (int j = 0; j < QK_K / 128; ++j) {
345 /* Load high bits */
346 const __m128i q4bitsH_0 = _mm_loadu_si128((const __m128i *)qh);
347 qh += 16;
348 const __m128i q4bitsH_1 = _mm_loadu_si128((const __m128i *)qh);
349 qh += 16;
350
351 /* Extract and shift high bits into position */
352 const __m128i q4h_0 = _mm_slli_epi16(_mm_and_si128(q4bitsH_0, m3), 4);
353 const __m128i q4h_1 = _mm_slli_epi16(_mm_and_si128(q4bitsH_1, m3), 4);
354 const __m128i q4h_2 = _mm_slli_epi16(_mm_and_si128(q4bitsH_0, _mm_set1_epi8(12)), 2);
355 const __m128i q4h_3 = _mm_slli_epi16(_mm_and_si128(q4bitsH_1, _mm_set1_epi8(12)), 2);
356 const __m128i q4h_4 = _mm_and_si128(q4bitsH_0, _mm_set1_epi8(48));
357 const __m128i q4h_5 = _mm_and_si128(q4bitsH_1, _mm_set1_epi8(48));
358 const __m128i q4h_6 = _mm_srli_epi16(_mm_and_si128(q4bitsH_0, _mm_set1_epi8(-64)), 2);
359 const __m128i q4h_7 = _mm_srli_epi16(_mm_and_si128(q4bitsH_1, _mm_set1_epi8(-64)), 2);
360
361 /* Load low bits */
362 const __m128i q4bits1_0 = _mm_loadu_si128((const __m128i *)ql);
363 ql += 16;
364 const __m128i q4bits1_1 = _mm_loadu_si128((const __m128i *)ql);
365 ql += 16;
366 const __m128i q4bits2_0 = _mm_loadu_si128((const __m128i *)ql);
367 ql += 16;
368 const __m128i q4bits2_1 = _mm_loadu_si128((const __m128i *)ql);
369 ql += 16;
370
371 /* Combine low and high bits to get 6-bit values (unsigned 0..63) */
372 const __m128i q4_0 = _mm_or_si128(_mm_and_si128(q4bits1_0, m15), q4h_0);
373 const __m128i q4_1 = _mm_or_si128(_mm_and_si128(q4bits1_1, m15), q4h_1);
374 const __m128i q4_2 = _mm_or_si128(_mm_and_si128(q4bits2_0, m15), q4h_2);
375 const __m128i q4_3 = _mm_or_si128(_mm_and_si128(q4bits2_1, m15), q4h_3);
376 const __m128i q4_4 = _mm_or_si128(_mm_and_si128(_mm_srli_epi16(q4bits1_0, 4), m15), q4h_4);
377 const __m128i q4_5 = _mm_or_si128(_mm_and_si128(_mm_srli_epi16(q4bits1_1, 4), m15), q4h_5);
378 const __m128i q4_6 = _mm_or_si128(_mm_and_si128(_mm_srli_epi16(q4bits2_0, 4), m15), q4h_6);
379 const __m128i q4_7 = _mm_or_si128(_mm_and_si128(_mm_srli_epi16(q4bits2_1, 4), m15), q4h_7);
380
381 /* Load Q8_K values */
382 const __m128i q8_0 = _mm_loadu_si128((const __m128i *)q8);
383 q8 += 16;
384 const __m128i q8_1 = _mm_loadu_si128((const __m128i *)q8);
385 q8 += 16;
386 const __m128i q8_2 = _mm_loadu_si128((const __m128i *)q8);
387 q8 += 16;
388 const __m128i q8_3 = _mm_loadu_si128((const __m128i *)q8);
389 q8 += 16;
390 const __m128i q8_4 = _mm_loadu_si128((const __m128i *)q8);
391 q8 += 16;
392 const __m128i q8_5 = _mm_loadu_si128((const __m128i *)q8);
393 q8 += 16;
394 const __m128i q8_6 = _mm_loadu_si128((const __m128i *)q8);
395 q8 += 16;
396 const __m128i q8_7 = _mm_loadu_si128((const __m128i *)q8);
397 q8 += 16;
398
399 /* Multiply: maddubs treats first arg as unsigned, second as signed */
400 __m128i p16_0 = _mm_maddubs_epi16(q4_0, q8_0);
401 __m128i p16_1 = _mm_maddubs_epi16(q4_1, q8_1);
402 __m128i p16_2 = _mm_maddubs_epi16(q4_2, q8_2);
403 __m128i p16_3 = _mm_maddubs_epi16(q4_3, q8_3);
404 __m128i p16_4 = _mm_maddubs_epi16(q4_4, q8_4);
405 __m128i p16_5 = _mm_maddubs_epi16(q4_5, q8_5);
406 __m128i p16_6 = _mm_maddubs_epi16(q4_6, q8_6);
407 __m128i p16_7 = _mm_maddubs_epi16(q4_7, q8_7);
408
409 /* Get scales for this iteration */
410 const __m128i scale_0 = _mm_shuffle_epi8(scales, _mm_loadu_si128((const __m128i *)q6k_scale_shuffle[is + 0]));
411 const __m128i scale_1 = _mm_shuffle_epi8(scales, _mm_loadu_si128((const __m128i *)q6k_scale_shuffle[is + 1]));
412 const __m128i scale_2 = _mm_shuffle_epi8(scales, _mm_loadu_si128((const __m128i *)q6k_scale_shuffle[is + 2]));
413 const __m128i scale_3 = _mm_shuffle_epi8(scales, _mm_loadu_si128((const __m128i *)q6k_scale_shuffle[is + 3]));
414 is += 4;
415
416 /* Scale the products and widen to 32-bit */
417 p16_0 = _mm_madd_epi16(_mm_cvtepi8_epi16(scale_0), p16_0);
418 p16_1 = _mm_madd_epi16(_mm_cvtepi8_epi16(_mm_bsrli_si128(scale_0, 8)), p16_1);
419 p16_2 = _mm_madd_epi16(_mm_cvtepi8_epi16(scale_1), p16_2);
420 p16_3 = _mm_madd_epi16(_mm_cvtepi8_epi16(_mm_bsrli_si128(scale_1, 8)), p16_3);
421 p16_4 = _mm_madd_epi16(_mm_cvtepi8_epi16(scale_2), p16_4);
422 p16_5 = _mm_madd_epi16(_mm_cvtepi8_epi16(_mm_bsrli_si128(scale_2, 8)), p16_5);
423 p16_6 = _mm_madd_epi16(_mm_cvtepi8_epi16(scale_3), p16_6);
424 p16_7 = _mm_madd_epi16(_mm_cvtepi8_epi16(_mm_bsrli_si128(scale_3, 8)), p16_7);
425
426 /* Accumulate */
427 sumi_0 = _mm_add_epi32(sumi_0, _mm_add_epi32(p16_0, p16_2));
428 sumi_1 = _mm_add_epi32(sumi_1, _mm_add_epi32(p16_1, p16_3));
429 sumi_0 = _mm_add_epi32(sumi_0, _mm_add_epi32(p16_4, p16_6));
430 sumi_1 = _mm_add_epi32(sumi_1, _mm_add_epi32(p16_5, p16_7));
431 }
432
433 /* Subtract the -32 offset contribution */
434 sumi_0 = _mm_sub_epi32(sumi_0, q8sclsub_0);
435 sumi_1 = _mm_sub_epi32(sumi_1, q8sclsub_1);
436
437 /* Combine and convert to float */
438 __m128i sumi = _mm_add_epi32(sumi_0, sumi_1);
439 __m128 sumf_vec = _mm_mul_ps(_mm_set1_ps(d), _mm_cvtepi32_ps(sumi));
440
441 /* Horizontal sum */
442 sumf_vec = _mm_hadd_ps(sumf_vec, sumf_vec);
443 sumf_vec = _mm_hadd_ps(sumf_vec, sumf_vec);
444 acc = _mm_add_ss(acc, sumf_vec);
445 }
446
447 return _mm_cvtss_f32(acc);
448}
449
450void gemv_q6_k_q8_k_sse(float *y,
451 const void *W,
452 const void *x_q8,
453 int M, int K)
454{
455 if (!y || !W || !x_q8 || M <= 0 || K <= 0) {
456 return;
457 }
458
459 const block_q6_K *blocks = (const block_q6_K *)W;
460 const block_q8_K *x = (const block_q8_K *)x_q8;
461 const int blocks_per_row = K / QK_K;
462
463 for (int row = 0; row < M; ++row) {
464 const block_q6_K *w_row = blocks + (size_t)row * (size_t)blocks_per_row;
465 y[row] = dot_q6_k_q8_k_sse(w_row, x, K);
466 }
467}
468#endif /* __SSSE3__ */
469
470/* ============================================================================
471 * AVX Implementation (for Sandy/Ivy Bridge - AVX without AVX2)
472 *
473 * Same as SSE but with prefetching for next block.
474 * Uses 128-bit integer ops (AVX doesn't add 256-bit int ops).
475 * ============================================================================ */
476
477#if defined(__AVX__) && !defined(__AVX2__)
478
479static float dot_q6_k_q8_k_avx(const block_q6_K *w,
480 const block_q8_K *x,
481 int K)
482{
483 const int nb = K / QK_K;
484 const __m128i m3 = _mm_set1_epi8(3);
485 const __m128i m15 = _mm_set1_epi8(15);
486
487 __m128 acc = _mm_setzero_ps();
488
489 for (int i = 0; i < nb; ++i) {
490 const float d = GGML_FP16_TO_FP32(w[i].d) * x[i].d;
491
492 const uint8_t *ql = w[i].ql;
493 const uint8_t *qh = w[i].qh;
494 const int8_t *q8 = x[i].qs;
495
496 /* Prefetch next block */
497 if (i + 1 < nb) {
498 _mm_prefetch((const char *)&w[i + 1], _MM_HINT_T0);
499 _mm_prefetch((const char *)&x[i + 1], _MM_HINT_T0);
500 }
501
502 /* Load scales and precompute the -32 offset contribution using bsums */
503 const __m128i scales = _mm_loadu_si128((const __m128i *)w[i].scales);
504 const __m128i q8sums_0 = _mm_loadu_si128((const __m128i *)x[i].bsums);
505 const __m128i q8sums_1 = _mm_loadu_si128((const __m128i *)x[i].bsums + 1);
506
507 /* Compute: sum(scale * bsum) * 32 for the -32 offset */
508 const __m128i scales_16_0 = _mm_cvtepi8_epi16(scales);
509 const __m128i scales_16_1 = _mm_cvtepi8_epi16(_mm_bsrli_si128(scales, 8));
510 const __m128i q8sclsub_0 = _mm_slli_epi32(_mm_madd_epi16(q8sums_0, scales_16_0), 5);
511 const __m128i q8sclsub_1 = _mm_slli_epi32(_mm_madd_epi16(q8sums_1, scales_16_1), 5);
512
513 __m128i sumi_0 = _mm_setzero_si128();
514 __m128i sumi_1 = _mm_setzero_si128();
515
516 int is = 0;
517
518 /* Process 256 weights in 2 iterations of 128 */
519 for (int j = 0; j < QK_K / 128; ++j) {
520 /* Load high bits */
521 const __m128i q4bitsH_0 = _mm_loadu_si128((const __m128i *)qh);
522 qh += 16;
523 const __m128i q4bitsH_1 = _mm_loadu_si128((const __m128i *)qh);
524 qh += 16;
525
526 /* Extract and shift high bits into position */
527 const __m128i q4h_0 = _mm_slli_epi16(_mm_and_si128(q4bitsH_0, m3), 4);
528 const __m128i q4h_1 = _mm_slli_epi16(_mm_and_si128(q4bitsH_1, m3), 4);
529 const __m128i q4h_2 = _mm_slli_epi16(_mm_and_si128(q4bitsH_0, _mm_set1_epi8(12)), 2);
530 const __m128i q4h_3 = _mm_slli_epi16(_mm_and_si128(q4bitsH_1, _mm_set1_epi8(12)), 2);
531 const __m128i q4h_4 = _mm_and_si128(q4bitsH_0, _mm_set1_epi8(48));
532 const __m128i q4h_5 = _mm_and_si128(q4bitsH_1, _mm_set1_epi8(48));
533 const __m128i q4h_6 = _mm_srli_epi16(_mm_and_si128(q4bitsH_0, _mm_set1_epi8(-64)), 2);
534 const __m128i q4h_7 = _mm_srli_epi16(_mm_and_si128(q4bitsH_1, _mm_set1_epi8(-64)), 2);
535
536 /* Load low bits */
537 const __m128i q4bits1_0 = _mm_loadu_si128((const __m128i *)ql);
538 ql += 16;
539 const __m128i q4bits1_1 = _mm_loadu_si128((const __m128i *)ql);
540 ql += 16;
541 const __m128i q4bits2_0 = _mm_loadu_si128((const __m128i *)ql);
542 ql += 16;
543 const __m128i q4bits2_1 = _mm_loadu_si128((const __m128i *)ql);
544 ql += 16;
545
546 /* Combine low and high bits to get 6-bit values (unsigned 0..63) */
547 const __m128i q4_0 = _mm_or_si128(_mm_and_si128(q4bits1_0, m15), q4h_0);
548 const __m128i q4_1 = _mm_or_si128(_mm_and_si128(q4bits1_1, m15), q4h_1);
549 const __m128i q4_2 = _mm_or_si128(_mm_and_si128(q4bits2_0, m15), q4h_2);
550 const __m128i q4_3 = _mm_or_si128(_mm_and_si128(q4bits2_1, m15), q4h_3);
551 const __m128i q4_4 = _mm_or_si128(_mm_and_si128(_mm_srli_epi16(q4bits1_0, 4), m15), q4h_4);
552 const __m128i q4_5 = _mm_or_si128(_mm_and_si128(_mm_srli_epi16(q4bits1_1, 4), m15), q4h_5);
553 const __m128i q4_6 = _mm_or_si128(_mm_and_si128(_mm_srli_epi16(q4bits2_0, 4), m15), q4h_6);
554 const __m128i q4_7 = _mm_or_si128(_mm_and_si128(_mm_srli_epi16(q4bits2_1, 4), m15), q4h_7);
555
556 /* Load Q8_K values */
557 const __m128i q8_0 = _mm_loadu_si128((const __m128i *)q8);
558 q8 += 16;
559 const __m128i q8_1 = _mm_loadu_si128((const __m128i *)q8);
560 q8 += 16;
561 const __m128i q8_2 = _mm_loadu_si128((const __m128i *)q8);
562 q8 += 16;
563 const __m128i q8_3 = _mm_loadu_si128((const __m128i *)q8);
564 q8 += 16;
565 const __m128i q8_4 = _mm_loadu_si128((const __m128i *)q8);
566 q8 += 16;
567 const __m128i q8_5 = _mm_loadu_si128((const __m128i *)q8);
568 q8 += 16;
569 const __m128i q8_6 = _mm_loadu_si128((const __m128i *)q8);
570 q8 += 16;
571 const __m128i q8_7 = _mm_loadu_si128((const __m128i *)q8);
572 q8 += 16;
573
574 /* Multiply: maddubs treats first arg as unsigned, second as signed */
575 __m128i p16_0 = _mm_maddubs_epi16(q4_0, q8_0);
576 __m128i p16_1 = _mm_maddubs_epi16(q4_1, q8_1);
577 __m128i p16_2 = _mm_maddubs_epi16(q4_2, q8_2);
578 __m128i p16_3 = _mm_maddubs_epi16(q4_3, q8_3);
579 __m128i p16_4 = _mm_maddubs_epi16(q4_4, q8_4);
580 __m128i p16_5 = _mm_maddubs_epi16(q4_5, q8_5);
581 __m128i p16_6 = _mm_maddubs_epi16(q4_6, q8_6);
582 __m128i p16_7 = _mm_maddubs_epi16(q4_7, q8_7);
583
584 /* Get scales for this iteration */
585 const __m128i scale_0 = _mm_shuffle_epi8(scales, _mm_loadu_si128((const __m128i *)q6k_scale_shuffle[is + 0]));
586 const __m128i scale_1 = _mm_shuffle_epi8(scales, _mm_loadu_si128((const __m128i *)q6k_scale_shuffle[is + 1]));
587 const __m128i scale_2 = _mm_shuffle_epi8(scales, _mm_loadu_si128((const __m128i *)q6k_scale_shuffle[is + 2]));
588 const __m128i scale_3 = _mm_shuffle_epi8(scales, _mm_loadu_si128((const __m128i *)q6k_scale_shuffle[is + 3]));
589 is += 4;
590
591 /* Scale the products and widen to 32-bit */
592 p16_0 = _mm_madd_epi16(_mm_cvtepi8_epi16(scale_0), p16_0);
593 p16_1 = _mm_madd_epi16(_mm_cvtepi8_epi16(_mm_bsrli_si128(scale_0, 8)), p16_1);
594 p16_2 = _mm_madd_epi16(_mm_cvtepi8_epi16(scale_1), p16_2);
595 p16_3 = _mm_madd_epi16(_mm_cvtepi8_epi16(_mm_bsrli_si128(scale_1, 8)), p16_3);
596 p16_4 = _mm_madd_epi16(_mm_cvtepi8_epi16(scale_2), p16_4);
597 p16_5 = _mm_madd_epi16(_mm_cvtepi8_epi16(_mm_bsrli_si128(scale_2, 8)), p16_5);
598 p16_6 = _mm_madd_epi16(_mm_cvtepi8_epi16(scale_3), p16_6);
599 p16_7 = _mm_madd_epi16(_mm_cvtepi8_epi16(_mm_bsrli_si128(scale_3, 8)), p16_7);
600
601 /* Accumulate */
602 sumi_0 = _mm_add_epi32(sumi_0, _mm_add_epi32(p16_0, p16_2));
603 sumi_1 = _mm_add_epi32(sumi_1, _mm_add_epi32(p16_1, p16_3));
604 sumi_0 = _mm_add_epi32(sumi_0, _mm_add_epi32(p16_4, p16_6));
605 sumi_1 = _mm_add_epi32(sumi_1, _mm_add_epi32(p16_5, p16_7));
606 }
607
608 /* Subtract the -32 offset contribution */
609 sumi_0 = _mm_sub_epi32(sumi_0, q8sclsub_0);
610 sumi_1 = _mm_sub_epi32(sumi_1, q8sclsub_1);
611
612 /* Combine and convert to float */
613 __m128i sumi = _mm_add_epi32(sumi_0, sumi_1);
614 __m128 sumf_vec = _mm_mul_ps(_mm_set1_ps(d), _mm_cvtepi32_ps(sumi));
615
616 /* Horizontal sum */
617 sumf_vec = _mm_hadd_ps(sumf_vec, sumf_vec);
618 sumf_vec = _mm_hadd_ps(sumf_vec, sumf_vec);
619 acc = _mm_add_ss(acc, sumf_vec);
620 }
621
622 return _mm_cvtss_f32(acc);
623}
624
625void gemv_q6_k_q8_k_avx(float *y,
626 const void *W,
627 const void *x_q8,
628 int M, int K)
629{
630 if (!y || !W || !x_q8 || M <= 0 || K <= 0) {
631 return;
632 }
633
634 const block_q6_K *blocks = (const block_q6_K *)W;
635 const block_q8_K *x = (const block_q8_K *)x_q8;
636 const int blocks_per_row = K / QK_K;
637
638 for (int row = 0; row < M; ++row) {
639 const block_q6_K *w_row = blocks + (size_t)row * (size_t)blocks_per_row;
640 y[row] = dot_q6_k_q8_k_avx(w_row, x, K);
641 }
642}
643
644#endif /* __AVX__ && !__AVX2__ */
645
646/* ============================================================================
647 * AVX2 Implementation (for modern CPUs with AVX2)
648 * ============================================================================ */
649
650#if defined(__AVX2__)
651
652/* Match ggml's x86 hsum_float_8 reduction tree exactly. Replacing the final
653 * two additions with _mm_hadd_ps changes one FP32 ULP for realistic MLP-down
654 * rows and can cross the absolute parity gate on AVX2/AVX-512 runners. */
655static inline float ck_q6k_hsum_float_8(const __m256 x)
656{
657 __m128 res = _mm256_extractf128_ps(x, 1);
658 res = _mm_add_ps(res, _mm256_castps256_ps128(x));
659 res = _mm_add_ps(res, _mm_movehl_ps(res, res));
660 res = _mm_add_ss(res, _mm_movehdup_ps(res));
661 return _mm_cvtss_f32(res);
662}
663
664/* Q6_K optimization note:
665 * ----------------------
666 * llama.cpp's newer Q6_K AVX2 work inspired an experiment that separates the
667 * zero-point correction from the hot unpack/dot loop by using Q8_K block sums
668 * (bsums). That can reduce repeated subtract-32 work and improve some CPU
669 * layouts, but it also changes reduction order and instruction balance.
670 *
671 * Current llama.cpp x86 kernels use the AVX2 reduction tree even in AVX-512
672 * builds. CK uses the same tree for its AVX2 and non-VBMI AVX-512 paths. Keep
673 * the final horizontal reduction in sync with ggml: an equivalent _mm_hadd_ps
674 * tree differs by one FP32 ULP on realistic MLP-down rows.
675 */
676
677static inline __m128i get_scale_shuffle_avx2(int i) {
678 static const uint8_t patterns[8][16] = {
679 { 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1 },
680 { 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3 },
681 { 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5 },
682 { 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7 },
683 { 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9 },
684 {10,10,10,10,10,10,10,10,11,11,11,11,11,11,11,11 },
685 {12,12,12,12,12,12,12,12,13,13,13,13,13,13,13,13 },
686 {14,14,14,14,14,14,14,14,15,15,15,15,15,15,15,15 },
687 };
688 return _mm_loadu_si128((const __m128i *)patterns[i]);
689}
690
691static float dot_q6_k_q8_k_avx2(const block_q6_K *w,
692 const block_q8_K *x,
693 int K)
694{
695 const int nb = K / QK_K;
696 const __m256i m3 = _mm256_set1_epi8(3);
697 const __m256i m15 = _mm256_set1_epi8(15);
698
699 __m256 acc = _mm256_setzero_ps();
700
701 for (int i = 0; i < nb; ++i) {
702 const float d = GGML_FP16_TO_FP32(w[i].d) * x[i].d;
703
704 const uint8_t *q4 = w[i].ql;
705 const uint8_t *qh = w[i].qh;
706 const int8_t *q8 = x[i].qs;
707
708 const __m256i q8sums = _mm256_loadu_si256((const __m256i *)x[i].bsums);
709 const __m128i scales = _mm_loadu_si128((const __m128i *)w[i].scales);
710 const __m256i scales_16 = _mm256_cvtepi8_epi16(scales);
711 const __m256i q8sclsub = _mm256_slli_epi32(
712 _mm256_madd_epi16(q8sums, scales_16), 5);
713
714 __m256i sumi = _mm256_setzero_si256();
715 int is = 0;
716
717 for (int j = 0; j < QK_K / 128; ++j) {
718 const __m256i q4bits1 = _mm256_loadu_si256((const __m256i *)q4);
719 q4 += 32;
720 const __m256i q4bits2 = _mm256_loadu_si256((const __m256i *)q4);
721 q4 += 32;
722 const __m256i q4bitsH = _mm256_loadu_si256((const __m256i *)qh);
723 qh += 32;
724
725 const __m256i q4h_0 = _mm256_slli_epi16(_mm256_and_si256(q4bitsH, m3), 4);
726 const __m256i q4h_1 = _mm256_slli_epi16(
727 _mm256_and_si256(q4bitsH, _mm256_set1_epi8(12)), 2);
728 const __m256i q4h_2 = _mm256_and_si256(q4bitsH, _mm256_set1_epi8(48));
729 const __m256i q4h_3 = _mm256_srli_epi16(
730 _mm256_and_si256(q4bitsH, _mm256_set1_epi8(-64)), 2);
731
732 const __m256i q4_0 = _mm256_or_si256(_mm256_and_si256(q4bits1, m15), q4h_0);
733 const __m256i q4_1 = _mm256_or_si256(_mm256_and_si256(q4bits2, m15), q4h_1);
734 const __m256i q4_2 = _mm256_or_si256(
735 _mm256_and_si256(_mm256_srli_epi16(q4bits1, 4), m15), q4h_2);
736 const __m256i q4_3 = _mm256_or_si256(
737 _mm256_and_si256(_mm256_srli_epi16(q4bits2, 4), m15), q4h_3);
738
739 const __m256i q8_0 = _mm256_loadu_si256((const __m256i *)q8);
740 q8 += 32;
741 const __m256i q8_1 = _mm256_loadu_si256((const __m256i *)q8);
742 q8 += 32;
743 const __m256i q8_2 = _mm256_loadu_si256((const __m256i *)q8);
744 q8 += 32;
745 const __m256i q8_3 = _mm256_loadu_si256((const __m256i *)q8);
746 q8 += 32;
747
748 __m256i p16_0 = _mm256_maddubs_epi16(q4_0, q8_0);
749 __m256i p16_1 = _mm256_maddubs_epi16(q4_1, q8_1);
750 __m256i p16_2 = _mm256_maddubs_epi16(q4_2, q8_2);
751 __m256i p16_3 = _mm256_maddubs_epi16(q4_3, q8_3);
752
753 const __m128i scale_0 = _mm_shuffle_epi8(
754 scales, get_scale_shuffle_avx2(is + 0));
755 const __m128i scale_1 = _mm_shuffle_epi8(
756 scales, get_scale_shuffle_avx2(is + 1));
757 const __m128i scale_2 = _mm_shuffle_epi8(
758 scales, get_scale_shuffle_avx2(is + 2));
759 const __m128i scale_3 = _mm_shuffle_epi8(
760 scales, get_scale_shuffle_avx2(is + 3));
761 is += 4;
762
763 p16_0 = _mm256_madd_epi16(_mm256_cvtepi8_epi16(scale_0), p16_0);
764 p16_1 = _mm256_madd_epi16(_mm256_cvtepi8_epi16(scale_1), p16_1);
765 p16_2 = _mm256_madd_epi16(_mm256_cvtepi8_epi16(scale_2), p16_2);
766 p16_3 = _mm256_madd_epi16(_mm256_cvtepi8_epi16(scale_3), p16_3);
767
768 sumi = _mm256_add_epi32(sumi, _mm256_add_epi32(p16_0, p16_1));
769 sumi = _mm256_add_epi32(sumi, _mm256_add_epi32(p16_2, p16_3));
770 }
771
772 sumi = _mm256_sub_epi32(sumi, q8sclsub);
773 acc = _mm256_fmadd_ps(_mm256_broadcast_ss(&d), _mm256_cvtepi32_ps(sumi), acc);
774 }
775
776 return ck_q6k_hsum_float_8(acc);
777}
778
779/*
780 * Four-token Q6_K x Q8_K microkernel.
781 *
782 * Each output keeps the same independent integer and FP32 reduction tree as
783 * dot_q6_k_q8_k_avx2(). The only changed ordering is between independent
784 * outputs: Q6 low/high-bit unpack and scale shuffles are performed once, then
785 * reused by as many as four Q8_K activation rows. This makes the provider
786 * eligible for exact llama.cpp parity while removing repeated Q6 decode work.
787 */
788static void dot_q6_k_q8_k_avx2_m4(const block_q6_K *w,
789 const block_q8_K *x,
790 int x_row_blocks,
791 int rows,
792 int K,
793 float out[4])
794{
795 const int nb = K / QK_K;
796 const __m256i m3 = _mm256_set1_epi8(3);
797 const __m256i m15 = _mm256_set1_epi8(15);
798 __m256 acc[4] = {
799 _mm256_setzero_ps(), _mm256_setzero_ps(),
800 _mm256_setzero_ps(), _mm256_setzero_ps()
801 };
802
803 for (int i = 0; i < nb; ++i) {
804 const uint8_t *q4 = w[i].ql;
805 const uint8_t *qh = w[i].qh;
806 const __m128i scales = _mm_loadu_si128((const __m128i *)w[i].scales);
807 const __m256i scales_16 = _mm256_cvtepi8_epi16(scales);
808 const float wd = GGML_FP16_TO_FP32(w[i].d);
809 __m256i sumi[4] = {
810 _mm256_setzero_si256(), _mm256_setzero_si256(),
811 _mm256_setzero_si256(), _mm256_setzero_si256()
812 };
813 int is = 0;
814
815 for (int j = 0; j < QK_K / 128; ++j) {
816 const __m256i q4bits1 = _mm256_loadu_si256((const __m256i *)q4);
817 q4 += 32;
818 const __m256i q4bits2 = _mm256_loadu_si256((const __m256i *)q4);
819 q4 += 32;
820 const __m256i q4bitsH = _mm256_loadu_si256((const __m256i *)qh);
821 qh += 32;
822
823 const __m256i q4h_0 = _mm256_slli_epi16(_mm256_and_si256(q4bitsH, m3), 4);
824 const __m256i q4h_1 = _mm256_slli_epi16(
825 _mm256_and_si256(q4bitsH, _mm256_set1_epi8(12)), 2);
826 const __m256i q4h_2 = _mm256_and_si256(q4bitsH, _mm256_set1_epi8(48));
827 const __m256i q4h_3 = _mm256_srli_epi16(
828 _mm256_and_si256(q4bitsH, _mm256_set1_epi8(-64)), 2);
829
830 const __m256i q6_0 = _mm256_or_si256(_mm256_and_si256(q4bits1, m15), q4h_0);
831 const __m256i q6_1 = _mm256_or_si256(_mm256_and_si256(q4bits2, m15), q4h_1);
832 const __m256i q6_2 = _mm256_or_si256(
833 _mm256_and_si256(_mm256_srli_epi16(q4bits1, 4), m15), q4h_2);
834 const __m256i q6_3 = _mm256_or_si256(
835 _mm256_and_si256(_mm256_srli_epi16(q4bits2, 4), m15), q4h_3);
836
837 const __m128i scale_0 = _mm_shuffle_epi8(
838 scales, get_scale_shuffle_avx2(is + 0));
839 const __m128i scale_1 = _mm_shuffle_epi8(
840 scales, get_scale_shuffle_avx2(is + 1));
841 const __m128i scale_2 = _mm_shuffle_epi8(
842 scales, get_scale_shuffle_avx2(is + 2));
843 const __m128i scale_3 = _mm_shuffle_epi8(
844 scales, get_scale_shuffle_avx2(is + 3));
845 is += 4;
846 const __m256i scale16_0 = _mm256_cvtepi8_epi16(scale_0);
847 const __m256i scale16_1 = _mm256_cvtepi8_epi16(scale_1);
848 const __m256i scale16_2 = _mm256_cvtepi8_epi16(scale_2);
849 const __m256i scale16_3 = _mm256_cvtepi8_epi16(scale_3);
850
851 for (int r = 0; r < rows; ++r) {
852 const block_q8_K *xb = x + (size_t)r * (size_t)x_row_blocks + i;
853 const int8_t *q8 = xb->qs + j * 128;
854 const __m256i q8_0 = _mm256_loadu_si256((const __m256i *)(q8 + 0));
855 const __m256i q8_1 = _mm256_loadu_si256((const __m256i *)(q8 + 32));
856 const __m256i q8_2 = _mm256_loadu_si256((const __m256i *)(q8 + 64));
857 const __m256i q8_3 = _mm256_loadu_si256((const __m256i *)(q8 + 96));
858
859 __m256i p16_0 = _mm256_maddubs_epi16(q6_0, q8_0);
860 __m256i p16_1 = _mm256_maddubs_epi16(q6_1, q8_1);
861 __m256i p16_2 = _mm256_maddubs_epi16(q6_2, q8_2);
862 __m256i p16_3 = _mm256_maddubs_epi16(q6_3, q8_3);
863 p16_0 = _mm256_madd_epi16(scale16_0, p16_0);
864 p16_1 = _mm256_madd_epi16(scale16_1, p16_1);
865 p16_2 = _mm256_madd_epi16(scale16_2, p16_2);
866 p16_3 = _mm256_madd_epi16(scale16_3, p16_3);
867
868 sumi[r] = _mm256_add_epi32(
869 sumi[r], _mm256_add_epi32(p16_0, p16_1));
870 sumi[r] = _mm256_add_epi32(
871 sumi[r], _mm256_add_epi32(p16_2, p16_3));
872 }
873 }
874
875 for (int r = 0; r < rows; ++r) {
876 const block_q8_K *xb = x + (size_t)r * (size_t)x_row_blocks + i;
877 const __m256i q8sums =
878 _mm256_loadu_si256((const __m256i *)xb->bsums);
879 const __m256i q8sclsub = _mm256_slli_epi32(
880 _mm256_madd_epi16(q8sums, scales_16), 5);
881 sumi[r] = _mm256_sub_epi32(sumi[r], q8sclsub);
882 const float d = wd * xb->d;
883 acc[r] = _mm256_fmadd_ps(
884 _mm256_broadcast_ss(&d), _mm256_cvtepi32_ps(sumi[r]), acc[r]);
885 }
886 }
887
888 for (int r = 0; r < rows; ++r) {
889 out[r] = ck_q6k_hsum_float_8(acc[r]);
890 }
891}
892
893static float dot_q6_k_prepared_q8_k_avx2(
894 const block_q6_K_prepared *w, const block_q8_K *x, int K)
895{
896 const int nb = K / QK_K;
897 __m256 acc = _mm256_setzero_ps();
898
899 for (int i = 0; i < nb; ++i) {
900 const __m128i scales =
901 _mm_loadu_si128((const __m128i *)(const void *)w[i].scales);
902 const __m256i scales_16 = _mm256_cvtepi8_epi16(scales);
903 const __m256i q8sums =
904 _mm256_loadu_si256((const __m256i *)(const void *)x[i].bsums);
905 const __m256i q8sclsub = _mm256_slli_epi32(
906 _mm256_madd_epi16(q8sums, scales_16), 5);
907 __m256i sumi = _mm256_setzero_si256();
908
909 for (int j = 0; j < 2; ++j) {
910 __m256i p16[4];
911 for (int lane = 0; lane < 4; ++lane) {
912 const int group = j * 4 + lane;
913 const __m256i q6 = _mm256_loadu_si256(
914 (const __m256i *)(const void *)(w[i].qs + group * 32));
915 const __m256i q8 = _mm256_loadu_si256(
916 (const __m256i *)(const void *)(x[i].qs + group * 32));
917 p16[lane] = _mm256_maddubs_epi16(q6, q8);
918 const __m128i scale = _mm_shuffle_epi8(
919 scales, get_scale_shuffle_avx2(group));
920 p16[lane] = _mm256_madd_epi16(
921 _mm256_cvtepi8_epi16(scale), p16[lane]);
922 }
923 sumi = _mm256_add_epi32(
924 sumi, _mm256_add_epi32(p16[0], p16[1]));
925 sumi = _mm256_add_epi32(
926 sumi, _mm256_add_epi32(p16[2], p16[3]));
927 }
928
929 sumi = _mm256_sub_epi32(sumi, q8sclsub);
930 const float d = GGML_FP16_TO_FP32(w[i].d) * x[i].d;
931 acc = _mm256_fmadd_ps(
932 _mm256_broadcast_ss(&d), _mm256_cvtepi32_ps(sumi), acc);
933 }
934
935 return ck_q6k_hsum_float_8(acc);
936}
937
938#if defined(__AVX512F__) && defined(__AVX512BW__) && \
939 defined(__AVX512VNNI__)
940/*
941 * Prepared Q6 stores unsigned 6-bit values as bytes, which maps directly to
942 * VPDPBUSD. Fold each pair of 32-byte groups into the same eight integer lanes
943 * used by the certified AVX2 provider. The correction, FP32 block FMA, and
944 * horizontal reduction therefore retain the established arithmetic contract.
945 */
946static float dot_q6_k_prepared_q8_k_avx512_vnni(
947 const block_q6_K_prepared *w, const block_q8_K *x, int K)
948{
949 const int nb = K / QK_K;
950 const __m512i scale_indices = _mm512_setr_epi32(
951 0, 0, 0, 0, 1, 1, 1, 1,
952 2, 2, 2, 2, 3, 3, 3, 3);
953 __m256 acc = _mm256_setzero_ps();
954
955 for (int i = 0; i < nb; ++i) {
956 const __m128i scales =
957 _mm_loadu_si128((const __m128i *)(const void *)w[i].scales);
958 const __m512i scales_32 = _mm512_cvtepi8_epi32(scales);
959 const __m256i scales_16 = _mm256_cvtepi8_epi16(scales);
960 const __m256i q8sums =
961 _mm256_loadu_si256((const __m256i *)(const void *)x[i].bsums);
962 const __m256i q8sclsub = _mm256_slli_epi32(
963 _mm256_madd_epi16(q8sums, scales_16), 5);
964 __m256i sumi = _mm256_setzero_si256();
965
966 for (int chunk = 0; chunk < 4; ++chunk) {
967 const __m512i q6 = _mm512_loadu_si512(
968 (const void *)(w[i].qs + chunk * 64));
969 const __m512i q8 = _mm512_loadu_si512(
970 (const void *)(x[i].qs + chunk * 64));
971 const __m512i products = _mm512_dpbusd_epi32(
972 _mm512_setzero_si512(), q6, q8);
973 const __m512i indices = _mm512_add_epi32(
974 scale_indices, _mm512_set1_epi32(chunk * 4));
975 const __m512i scaled = _mm512_mullo_epi32(
976 products, _mm512_permutexvar_epi32(indices, scales_32));
977 const __m512i upper_half = _mm512_shuffle_i32x4(
978 scaled, scaled, _MM_SHUFFLE(3, 2, 3, 2));
979 const __m256i pair_sum = _mm256_add_epi32(
980 _mm512_castsi512_si256(scaled),
981 _mm512_castsi512_si256(upper_half));
982 sumi = _mm256_add_epi32(sumi, pair_sum);
983 }
984
985 sumi = _mm256_sub_epi32(sumi, q8sclsub);
986 const float d = GGML_FP16_TO_FP32(w[i].d) * x[i].d;
987 acc = _mm256_fmadd_ps(
988 _mm256_broadcast_ss(&d), _mm256_cvtepi32_ps(sumi), acc);
989 }
990
991 return ck_q6k_hsum_float_8(acc);
992}
993#endif
994
995void gemv_q6_k_q8_k_avx2(float *y,
996 const void *W,
997 const void *x_q8,
998 int M, int K)
999{
1000 if (!y || !W || !x_q8 || M <= 0 || K <= 0) {
1001 return;
1002 }
1003
1004 const block_q6_K *blocks = (const block_q6_K *)W;
1005 const block_q8_K *x = (const block_q8_K *)x_q8;
1006 const int blocks_per_row = K / QK_K;
1007
1008 for (int row = 0; row < M; ++row) {
1009 const block_q6_K *w_row = blocks + (size_t)row * (size_t)blocks_per_row;
1010 y[row] = dot_q6_k_q8_k_avx2(w_row, x, K);
1011 }
1012}
1013#endif /* __AVX2__ */
1014
1015/* ============================================================================
1016 * AVX-512 Implementation
1017 *
1018 * Uses 512-bit ZMM registers to process 64 bytes at a time.
1019 * Processes entire 256-element Q6_K block in fewer iterations.
1020 * ============================================================================ */
1021
1022#if defined(__AVX512F__) && defined(__AVX512BW__) && defined(__AVX512VBMI__)
1023
1024/**
1025 * @brief AVX-512 dot product for Q6_K x Q8_K with VBMI
1026 *
1027 * Uses AVX-512 VBMI for efficient byte permutation.
1028 */
1029static float dot_q6_k_q8_k_avx512_vbmi(const block_q6_K *w,
1030 const block_q8_K *x,
1031 int K)
1032{
1033 const int nb = K / QK_K;
1034 const __m512i m4 = _mm512_set1_epi8(0xF);
1035 const __m512i m2 = _mm512_set1_epi8(3);
1036 const __m512i m32s = _mm512_set1_epi8(32);
1037
1038 __m512 acc = _mm512_setzero_ps();
1039
1040 for (int i = 0; i < nb; ++i) {
1041 const float d = GGML_FP16_TO_FP32(w[i].d) * x[i].d;
1042
1043 const uint8_t *ql = w[i].ql;
1044 const uint8_t *qh = w[i].qh;
1045 const int8_t *q8 = x[i].qs;
1046 const int8_t *sc = w[i].scales;
1047
1048 __m512i sumi = _mm512_setzero_si512();
1049
1050 /* Process 256 weights in one iteration using AVX-512 */
1051 /* Load 64 bytes of low bits (covers 128 weights, need 2 loads for full block) */
1052 const __m512i q4bits1 = _mm512_loadu_si512((const __m512i *)ql); /* ql[0..63] */
1053 const __m512i q4bits2 = _mm512_loadu_si512((const __m512i *)(ql + 64)); /* ql[64..127] */
1054
1055 /* Load 64 bytes of high bits */
1056 const __m512i q4bitsH = _mm512_loadu_si512((const __m512i *)qh);
1057
1058 /* Extract high 2-bit contributions for each group of 32 weights */
1059 /* Group 0: bits 0-1 of qh -> weights 0-31 */
1060 const __m512i q4h_0 = _mm512_slli_epi16(_mm512_and_si512(q4bitsH, m2), 4);
1061 /* Group 1: bits 2-3 of qh -> weights 32-63 */
1062 const __m512i q4h_1 = _mm512_slli_epi16(_mm512_and_si512(_mm512_srli_epi16(q4bitsH, 2), m2), 4);
1063 /* Group 2: bits 4-5 of qh -> weights 64-95 */
1064 const __m512i q4h_2 = _mm512_slli_epi16(_mm512_and_si512(_mm512_srli_epi16(q4bitsH, 4), m2), 4);
1065 /* Group 3: bits 6-7 of qh -> weights 96-127 */
1066 const __m512i q4h_3 = _mm512_slli_epi16(_mm512_and_si512(_mm512_srli_epi16(q4bitsH, 6), m2), 4);
1067
1068 /* Combine low nibbles with high bits to get 6-bit values (0-63) */
1069 /* First 64 weights: low nibbles of ql[0..63] */
1070 const __m512i q6_0 = _mm512_or_si512(_mm512_and_si512(q4bits1, m4), q4h_0);
1071 const __m512i q6_1 = _mm512_or_si512(_mm512_and_si512(q4bits2, m4), q4h_1);
1072 /* Second 64 weights: high nibbles of ql[0..63] */
1073 const __m512i q6_2 = _mm512_or_si512(_mm512_and_si512(_mm512_srli_epi16(q4bits1, 4), m4), q4h_2);
1074 const __m512i q6_3 = _mm512_or_si512(_mm512_and_si512(_mm512_srli_epi16(q4bits2, 4), m4), q4h_3);
1075
1076 /* Load Q8_K values (256 int8 values = 4 x 64) */
1077 const __m512i q8_0 = _mm512_loadu_si512((const __m512i *)q8);
1078 const __m512i q8_1 = _mm512_loadu_si512((const __m512i *)(q8 + 64));
1079 const __m512i q8_2 = _mm512_loadu_si512((const __m512i *)(q8 + 128));
1080 const __m512i q8_3 = _mm512_loadu_si512((const __m512i *)(q8 + 192));
1081
1082 /* Compute 32 * q8 for the offset subtraction */
1083 __m512i q8s_0 = _mm512_maddubs_epi16(m32s, q8_0);
1084 __m512i q8s_1 = _mm512_maddubs_epi16(m32s, q8_1);
1085 __m512i q8s_2 = _mm512_maddubs_epi16(m32s, q8_2);
1086 __m512i q8s_3 = _mm512_maddubs_epi16(m32s, q8_3);
1087
1088 /* Multiply unsigned q6 * signed q8 */
1089 __m512i p16_0 = _mm512_maddubs_epi16(q6_0, q8_0);
1090 __m512i p16_1 = _mm512_maddubs_epi16(q6_1, q8_1);
1091 __m512i p16_2 = _mm512_maddubs_epi16(q6_2, q8_2);
1092 __m512i p16_3 = _mm512_maddubs_epi16(q6_3, q8_3);
1093
1094 /* Subtract offset: (q6 - 32) * q8 = q6*q8 - 32*q8 */
1095 p16_0 = _mm512_sub_epi16(p16_0, q8s_0);
1096 p16_1 = _mm512_sub_epi16(p16_1, q8s_1);
1097 p16_2 = _mm512_sub_epi16(p16_2, q8s_2);
1098 p16_3 = _mm512_sub_epi16(p16_3, q8s_3);
1099
1100 /* Load and broadcast scales using VBMI permute
1101 * Each scale applies to 16 weights, so we need to broadcast appropriately
1102 * scales[0..15] for the 16 sub-blocks */
1103 const __m128i scales_128 = _mm_loadu_si128((const __m128i *)sc);
1104
1105 /* Create scale broadcast patterns for 64 weights (4 scales per 64 weights) */
1106 /* Pattern: each scale repeated 16 times for 16 weights */
1107 const __m512i scale_idx_0 = _mm512_set_epi8(
1108 3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
1109 2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
1110 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1111 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
1112 const __m512i scale_idx_1 = _mm512_set_epi8(
1113 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
1114 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
1115 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
1116 4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4);
1117 const __m512i scale_idx_2 = _mm512_set_epi8(
1118 11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,
1119 10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,
1120 9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,
1121 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8);
1122 const __m512i scale_idx_3 = _mm512_set_epi8(
1123 15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,
1124 14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,
1125 13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,
1126 12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12);
1127
1128 /* Broadcast scales to 512-bit using VBMI permutexvar */
1129 const __m512i scales_512 = _mm512_broadcast_i32x4(scales_128);
1130 const __m512i sc_0 = _mm512_permutexvar_epi8(scale_idx_0, scales_512);
1131 const __m512i sc_1 = _mm512_permutexvar_epi8(scale_idx_1, scales_512);
1132 const __m512i sc_2 = _mm512_permutexvar_epi8(scale_idx_2, scales_512);
1133 const __m512i sc_3 = _mm512_permutexvar_epi8(scale_idx_3, scales_512);
1134
1135 /* Sign-extend scales to 16-bit and multiply with products */
1136 /* For efficiency, we process in two halves (low and high 256 bits) */
1137 __m512i p32_0 = _mm512_madd_epi16(_mm512_cvtepi8_epi16(_mm512_castsi512_si256(sc_0)), p16_0);
1138 __m512i p32_1 = _mm512_madd_epi16(_mm512_cvtepi8_epi16(_mm512_castsi512_si256(sc_1)), p16_1);
1139 __m512i p32_2 = _mm512_madd_epi16(_mm512_cvtepi8_epi16(_mm512_castsi512_si256(sc_2)), p16_2);
1140 __m512i p32_3 = _mm512_madd_epi16(_mm512_cvtepi8_epi16(_mm512_castsi512_si256(sc_3)), p16_3);
1141
1142 /* Accumulate */
1143 sumi = _mm512_add_epi32(sumi, p32_0);
1144 sumi = _mm512_add_epi32(sumi, p32_1);
1145 sumi = _mm512_add_epi32(sumi, p32_2);
1146 sumi = _mm512_add_epi32(sumi, p32_3);
1147
1148 /* Scale by d and accumulate */
1149 acc = _mm512_fmadd_ps(_mm512_set1_ps(d), _mm512_cvtepi32_ps(sumi), acc);
1150 }
1151
1152 return _mm512_reduce_add_ps(acc);
1153}
1154
1155void gemv_q6_k_q8_k_avx512_vbmi(float *y,
1156 const void *W,
1157 const void *x_q8,
1158 int M, int K)
1159{
1160 if (!y || !W || !x_q8 || M <= 0 || K <= 0) {
1161 return;
1162 }
1163
1164 const block_q6_K *blocks = (const block_q6_K *)W;
1165 const block_q8_K *x = (const block_q8_K *)x_q8;
1166 const int blocks_per_row = K / QK_K;
1167
1168 for (int row = 0; row < M; ++row) {
1169 const block_q6_K *w_row = blocks + (size_t)row * (size_t)blocks_per_row;
1170 y[row] = dot_q6_k_q8_k_avx512_vbmi(w_row, x, K);
1171 }
1172}
1173
1174#endif /* __AVX512F__ && __AVX512BW__ && __AVX512VBMI__ */
1175
1176#if defined(__AVX512F__) && defined(__AVX512BW__)
1177
1178/**
1179 * @brief AVX-512 dot product for Q6_K x Q8_K
1180 *
1181 * Works on all AVX-512 CPUs (Skylake-X and newer).
1182 * Uses same algorithm as AVX2, but benefits from AVX-512's wider FMA
1183 * and efficient horizontal reduction.
1184 */
1185static float dot_q6_k_q8_k_avx512(const block_q6_K *w,
1186 const block_q8_K *x,
1187 int K)
1188{
1189 const int nb = K / QK_K;
1190 const __m256i m4 = _mm256_set1_epi8(0xF);
1191 const __m256i m2 = _mm256_set1_epi8(3);
1192 const __m256i m32s = _mm256_set1_epi8(32);
1193
1194 /* Use 256-bit float accumulator, same as AVX2 */
1195 __m256 acc = _mm256_setzero_ps();
1196
1197 for (int i = 0; i < nb; ++i) {
1198 const float d = GGML_FP16_TO_FP32(w[i].d) * x[i].d;
1199
1200 const uint8_t *q4 = w[i].ql;
1201 const uint8_t *qh = w[i].qh;
1202 const int8_t *q8 = x[i].qs;
1203
1204 const __m128i scales = _mm_loadu_si128((const __m128i *)w[i].scales);
1205
1206 /* Use 256-bit integer accumulator, same as AVX2 */
1207 __m256i sumi = _mm256_setzero_si256();
1208 int is = 0;
1209
1210 /* Process 256 weights in 2 iterations of 128 (same structure as AVX2) */
1211 for (int j = 0; j < QK_K / 128; ++j) {
1212 /* Get scale shuffle patterns - identical to AVX2 */
1213 static const uint8_t patterns[8][16] = {
1214 { 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1 },
1215 { 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3 },
1216 { 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5 },
1217 { 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7 },
1218 { 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9 },
1219 {10,10,10,10,10,10,10,10,11,11,11,11,11,11,11,11 },
1220 {12,12,12,12,12,12,12,12,13,13,13,13,13,13,13,13 },
1221 {14,14,14,14,14,14,14,14,15,15,15,15,15,15,15,15 },
1222 };
1223
1224 const __m128i scale_0 = _mm_shuffle_epi8(scales, _mm_loadu_si128((const __m128i *)patterns[is + 0]));
1225 const __m128i scale_1 = _mm_shuffle_epi8(scales, _mm_loadu_si128((const __m128i *)patterns[is + 1]));
1226 const __m128i scale_2 = _mm_shuffle_epi8(scales, _mm_loadu_si128((const __m128i *)patterns[is + 2]));
1227 const __m128i scale_3 = _mm_shuffle_epi8(scales, _mm_loadu_si128((const __m128i *)patterns[is + 3]));
1228 is += 4;
1229
1230 /* Load low bits */
1231 const __m256i q4bits1 = _mm256_loadu_si256((const __m256i *)q4);
1232 q4 += 32;
1233 const __m256i q4bits2 = _mm256_loadu_si256((const __m256i *)q4);
1234 q4 += 32;
1235 const __m256i q4bitsH = _mm256_loadu_si256((const __m256i *)qh);
1236 qh += 32;
1237
1238 /* Extract high 2-bit contributions */
1239 const __m256i q4h_0 = _mm256_slli_epi16(_mm256_and_si256(q4bitsH, m2), 4);
1240 const __m256i q4h_1 = _mm256_slli_epi16(_mm256_and_si256(_mm256_srli_epi16(q4bitsH, 2), m2), 4);
1241 const __m256i q4h_2 = _mm256_slli_epi16(_mm256_and_si256(_mm256_srli_epi16(q4bitsH, 4), m2), 4);
1242 const __m256i q4h_3 = _mm256_slli_epi16(_mm256_and_si256(_mm256_srli_epi16(q4bitsH, 6), m2), 4);
1243
1244 /* Combine low + high bits to get 6-bit values */
1245 const __m256i q4_0 = _mm256_or_si256(_mm256_and_si256(q4bits1, m4), q4h_0);
1246 const __m256i q4_1 = _mm256_or_si256(_mm256_and_si256(q4bits2, m4), q4h_1);
1247 const __m256i q4_2 = _mm256_or_si256(_mm256_and_si256(_mm256_srli_epi16(q4bits1, 4), m4), q4h_2);
1248 const __m256i q4_3 = _mm256_or_si256(_mm256_and_si256(_mm256_srli_epi16(q4bits2, 4), m4), q4h_3);
1249
1250 /* Load Q8_K values */
1251 const __m256i q8_0 = _mm256_loadu_si256((const __m256i *)q8);
1252 q8 += 32;
1253 const __m256i q8_1 = _mm256_loadu_si256((const __m256i *)q8);
1254 q8 += 32;
1255 const __m256i q8_2 = _mm256_loadu_si256((const __m256i *)q8);
1256 q8 += 32;
1257 const __m256i q8_3 = _mm256_loadu_si256((const __m256i *)q8);
1258 q8 += 32;
1259
1260 /* Compute 32 * q8 for offset */
1261 __m256i q8s_0 = _mm256_maddubs_epi16(m32s, q8_0);
1262 __m256i q8s_1 = _mm256_maddubs_epi16(m32s, q8_1);
1263 __m256i q8s_2 = _mm256_maddubs_epi16(m32s, q8_2);
1264 __m256i q8s_3 = _mm256_maddubs_epi16(m32s, q8_3);
1265
1266 /* Multiply q4 * q8 (unsigned * signed) */
1267 __m256i p16_0 = _mm256_maddubs_epi16(q4_0, q8_0);
1268 __m256i p16_1 = _mm256_maddubs_epi16(q4_1, q8_1);
1269 __m256i p16_2 = _mm256_maddubs_epi16(q4_2, q8_2);
1270 __m256i p16_3 = _mm256_maddubs_epi16(q4_3, q8_3);
1271
1272 /* Subtract offset: (q4 - 32) * q8 = q4*q8 - 32*q8 */
1273 p16_0 = _mm256_sub_epi16(p16_0, q8s_0);
1274 p16_1 = _mm256_sub_epi16(p16_1, q8s_1);
1275 p16_2 = _mm256_sub_epi16(p16_2, q8s_2);
1276 p16_3 = _mm256_sub_epi16(p16_3, q8s_3);
1277
1278 /* Apply scales - produces 8 int32 each */
1279 p16_0 = _mm256_madd_epi16(_mm256_cvtepi8_epi16(scale_0), p16_0);
1280 p16_1 = _mm256_madd_epi16(_mm256_cvtepi8_epi16(scale_1), p16_1);
1281 p16_2 = _mm256_madd_epi16(_mm256_cvtepi8_epi16(scale_2), p16_2);
1282 p16_3 = _mm256_madd_epi16(_mm256_cvtepi8_epi16(scale_3), p16_3);
1283
1284 /* Accumulate all 4 into sumi (same as AVX2) */
1285 sumi = _mm256_add_epi32(sumi, _mm256_add_epi32(p16_0, p16_1));
1286 sumi = _mm256_add_epi32(sumi, _mm256_add_epi32(p16_2, p16_3));
1287 }
1288
1289 /* Scale by d and accumulate */
1290 acc = _mm256_fmadd_ps(_mm256_set1_ps(d), _mm256_cvtepi32_ps(sumi), acc);
1291 }
1292
1293 return ck_q6k_hsum_float_8(acc);
1294}
1295
1296void gemv_q6_k_q8_k_avx512(float *y,
1297 const void *W,
1298 const void *x_q8,
1299 int M, int K)
1300{
1301 if (!y || !W || !x_q8 || M <= 0 || K <= 0) {
1302 return;
1303 }
1304
1305 const block_q6_K *blocks = (const block_q6_K *)W;
1306 const block_q8_K *x = (const block_q8_K *)x_q8;
1307 const int blocks_per_row = K / QK_K;
1308
1309 for (int row = 0; row < M; ++row) {
1310 const block_q6_K *w_row = blocks + (size_t)row * (size_t)blocks_per_row;
1311 y[row] = dot_q6_k_q8_k_avx512(w_row, x, K);
1312 }
1313}
1314
1315#endif /* __AVX512F__ && __AVX512BW__ */
1316
1317/* ============================================================================
1318 * Dispatch Functions
1319 * ============================================================================ */
1320
1321/**
1322 * @brief Q6_K x Q8_K dot product (single row)
1323 */
1324void vec_dot_q6_k_q8_k(int n, float *s, const void *vx, const void *vy)
1325{
1326 if (!s || !vx || !vy || n <= 0) {
1327 return;
1328 }
1329
1330 const block_q6_K *x = (const block_q6_K *)vx;
1331 const block_q8_K *y = (const block_q8_K *)vy;
1332
1333 /* This is the architecture-neutral scalar oracle. Production x86 dispatch
1334 * uses the separately parity-tested SIMD reduction tree. */
1335 *s = dot_q6_k_q8_k_ref(x, y, n);
1336}
1337
1338/**
1339 * @brief GEMV: y = W @ x where W is Q6_K and x is Q8_K
1340 */
1341void gemv_q6_k_q8_k(float *y,
1342 const void *W,
1343 const void *x_q8,
1344 int M, int K)
1345{
1347 gemv_q6_k_q8_k_ref(y, W, x_q8, M, K);
1348 return;
1349 }
1350
1351#if defined(__AVX2__)
1352 /* llama.cpp's x86 Q6_K production graph keeps the AVX2 reduction order
1353 * even when AVX-512 is available. Wider ISA availability is not a license
1354 * to change this numerical contract. */
1355 gemv_q6_k_q8_k_avx2(y, W, x_q8, M, K);
1356 return;
1357#elif defined(__AVX__)
1358 gemv_q6_k_q8_k_avx(y, W, x_q8, M, K);
1359 return;
1360#elif defined(__SSE4_1__)
1361 gemv_q6_k_q8_k_sse(y, W, x_q8, M, K);
1362 return;
1363#endif
1364 gemv_q6_k_q8_k_ref(y, W, x_q8, M, K);
1365}
1366
1368{
1370 return "q6_k_q8_k_ref";
1371 }
1372#if defined(__AVX2__)
1373 return "q6_k_q8_k_avx2";
1374#elif defined(__AVX__)
1375 return "q6_k_q8_k_avx";
1376#elif defined(__SSE4_1__)
1377 return "q6_k_q8_k_sse";
1378#else
1379 return "q6_k_q8_k_ref";
1380#endif
1381}
1382
1383/* ============================================================================
1384 * PARALLEL VERSIONS (for parallel orchestration)
1385 *
1386 * These receive ith (thread index) and nth (total threads) from orchestration.
1387 * OpenMP lives in orchestration layer, NOT here.
1388 *
1389 * Naming: *_parallel = receives ith/nth, processes only its portion
1390 * ============================================================================ */
1391
1392/**
1393 * @brief Parallel reference GEMV for Q6_K × Q8_K
1394 *
1395 * Caller provides ith (thread index) and nth (total threads).
1396 * Each thread processes rows [r0, r1).
1397 */
1399 const void *W,
1400 const void *x_q8,
1401 int M, int K,
1402 int ith, int nth)
1403{
1404 if (!y || !W || !x_q8 || M <= 0 || K <= 0) return;
1405 if (ith < 0 || nth <= 0 || ith >= nth) return;
1406
1407 /* Compute row range for this thread */
1408 const int dr = (M + nth - 1) / nth;
1409 const int r0 = dr * ith;
1410 const int r1 = (r0 + dr < M) ? (r0 + dr) : M;
1411
1412 if (r0 >= M) return;
1413
1414 const block_q6_K *blocks = (const block_q6_K *)W;
1415 const block_q8_K *x = (const block_q8_K *)x_q8;
1416 const int blocks_per_row = K / QK_K;
1417
1418 for (int row = r0; row < r1; ++row) {
1419 const block_q6_K *w_row = blocks + (size_t)row * (size_t)blocks_per_row;
1420 y[row] = dot_q6_k_q8_k_ref(w_row, x, K);
1421 }
1422}
1423
1424/**
1425 * @brief Parallel SIMD GEMV for Q6_K × Q8_K
1426 *
1427 * Uses best available SIMD (AVX/SSE) with row prefetching.
1428 * Caller provides ith/nth from OpenMP region.
1429 */
1431 const void *W,
1432 const void *x_q8,
1433 int M, int K,
1434 int ith, int nth)
1435{
1436 if (!y || !W || !x_q8 || M <= 0 || K <= 0) return;
1437 if (ith < 0 || nth <= 0 || ith >= nth) return;
1438
1439 const int dr = (M + nth - 1) / nth;
1440 const int r0 = dr * ith;
1441 const int r1 = (r0 + dr < M) ? (r0 + dr) : M;
1442
1443 if (r0 >= M) return;
1444
1445 const block_q6_K *blocks = (const block_q6_K *)W;
1446 const block_q8_K *x = (const block_q8_K *)x_q8;
1447 const int blocks_per_row = K / QK_K;
1448 const int strict = ck_strict_parity_enabled() || ck_q6k_q8k_force_ref();
1449
1450 for (int row = r0; row < r1; ++row) {
1451 const block_q6_K *w_row = blocks + (size_t)row * (size_t)blocks_per_row;
1452#if defined(__AVX2__)
1453 y[row] = strict ? dot_q6_k_q8_k_ref(w_row, x, K)
1454 : dot_q6_k_q8_k_avx2(w_row, x, K);
1455#elif defined(__AVX__)
1456 y[row] = strict ? dot_q6_k_q8_k_ref(w_row, x, K)
1457 : dot_q6_k_q8_k_avx(w_row, x, K);
1458#elif defined(__SSE4_1__)
1459 y[row] = strict ? dot_q6_k_q8_k_ref(w_row, x, K)
1460 : dot_q6_k_q8_k_sse(w_row, x, K);
1461#else
1462 y[row] = dot_q6_k_q8_k_ref(w_row, x, K);
1463#endif
1464 }
1465}
1466
1467static inline float ck_dot_q6_k_q8_k_fast_or_ref(const block_q6_K *w,
1468 const block_q8_K *x,
1469 int K);
1470
1471/**
1472 * @brief GEMM: Y = W @ X^T where W is Q6_K and X is Q8_K
1473 *
1474 * @param Y Output matrix [N x M] in row-major
1475 * @param W Weight matrix in Q6_K format [M x K]
1476 * @param X_q8 Input matrix in Q8_K format [N x K]
1477 * @param M Number of output rows (output dim)
1478 * @param N Number of input vectors (batch size)
1479 * @param K Input dimension
1480 */
1481void gemm_q6_k_q8_k(float *Y,
1482 const void *W,
1483 const void *X_q8,
1484 int M, int N, int K)
1485{
1486 if (!Y || !W || !X_q8 || M <= 0 || N <= 0 || K <= 0) {
1487 return;
1488 }
1489
1490 const block_q8_K *X = (const block_q8_K *)X_q8;
1491 const int blocks_per_vec = K / QK_K;
1492
1493 for (int n = 0; n < N; ++n) {
1494 const block_q8_K *x_row = X + (size_t)n * (size_t)blocks_per_vec;
1495 gemv_q6_k_q8_k(&Y[n * M], W, x_row, M, K);
1496 }
1497}
1498
1499/**
1500 * @brief NT GEMM: C = A @ B^T where A is Q8_K and B is Q6_K
1501 *
1502 * This is the typical inference pattern:
1503 * - A: Activations in Q8_K format [M x K]
1504 * - B: Weights in Q6_K format [N x K]
1505 * - C: Output [M x N]
1506 *
1507 * @param A_q8 Input activations in Q8_K format
1508 * @param B Weight matrix in Q6_K format
1509 * @param bias Optional bias vector [N]
1510 * @param C Output matrix
1511 * @param M Batch size (number of tokens)
1512 * @param N Output dimension
1513 * @param K Input dimension
1514 */
1515void gemm_nt_q6_k_q8_k(const void *A_q8,
1516 const void *B,
1517 const float *bias,
1518 float *C,
1519 int M, int N, int K)
1520{
1521 if (!A_q8 || !B || !C) {
1522 return;
1523 }
1524 if (M <= 0 || N <= 0 || K <= 0) {
1525 return;
1526 }
1527
1528 /* Prefill GEMM is the hot Qwen2/Qwen3.5 MLP-down path. Keep decode
1529 * gemv_q6_k_q8_k() conservative, but allow GEMM/prefill to use the
1530 * parity-gated SIMD dot helper by default. CK strict parity and
1531 * CK_DEBUG_Q6K_Q8K_REF=1 still force the scalar reference reduction. */
1532 const block_q8_K *A = (const block_q8_K *)A_q8;
1533 const block_q6_K *W = (const block_q6_K *)B;
1534 const int blocks_per_vec = K / QK_K;
1535 const int blocks_per_row = K / QK_K;
1536
1537 for (int m = 0; m < M; ++m) {
1538 const block_q8_K *a_row = A + (size_t)m * (size_t)blocks_per_vec;
1539 float *c_row = C + (size_t)m * (size_t)N;
1540 for (int n = 0; n < N; ++n) {
1541 const block_q6_K *w_row = W + (size_t)n * (size_t)blocks_per_row;
1542 const float b = bias ? bias[n] : 0.0f;
1543 c_row[n] = ck_dot_q6_k_q8_k_fast_or_ref(w_row, a_row, K) + b;
1544 }
1545 }
1546}
1547
1549 const void *A_q8,
1550 const void *B_prepared,
1551 const float *bias,
1552 float *C,
1553 int M, int N, int K,
1554 int m0, int m1,
1555 int n0, int n1,
1556 int use_avx512_vnni)
1557{
1558#if !defined(__AVX2__)
1559 (void)A_q8; (void)B_prepared; (void)bias; (void)C;
1560 (void)M; (void)N; (void)K; (void)m0; (void)m1; (void)n0; (void)n1;
1561 (void)use_avx512_vnni;
1562#else
1563#if !defined(__AVX512F__) || !defined(__AVX512BW__) || \
1564 !defined(__AVX512VNNI__)
1565 (void)use_avx512_vnni;
1566#endif
1567 if (!A_q8 || !B_prepared || !C || M <= 0 || N <= 0 || K <= 0 ||
1568 (K % QK_K) != 0) return;
1569 if (m0 < 0) m0 = 0;
1570 if (n0 < 0) n0 = 0;
1571 if (m1 > M) m1 = M;
1572 if (n1 > N) n1 = N;
1573 if (m0 >= m1 || n0 >= n1) return;
1574
1575 const block_q8_K *A = (const block_q8_K *)A_q8;
1576 const block_q6_K_prepared *W =
1577 (const block_q6_K_prepared *)B_prepared;
1578 const int blocks_per_row = K / QK_K;
1579 for (int n = n0; n < n1; ++n) {
1580 const block_q6_K_prepared *w_row =
1581 W + (size_t)n * (size_t)blocks_per_row;
1582 const float b = bias ? bias[n] : 0.0f;
1583 for (int m = m0; m < m1; ++m) {
1584 const block_q8_K *a_row =
1585 A + (size_t)m * (size_t)blocks_per_row;
1586#if defined(__AVX512F__) && defined(__AVX512BW__) && \
1587 defined(__AVX512VNNI__)
1588 const float dot = use_avx512_vnni
1589 ? dot_q6_k_prepared_q8_k_avx512_vnni(w_row, a_row, K)
1590 : dot_q6_k_prepared_q8_k_avx2(w_row, a_row, K);
1591#else
1592 const float dot = dot_q6_k_prepared_q8_k_avx2(w_row, a_row, K);
1593#endif
1594 C[(size_t)m * (size_t)N + (size_t)n] =
1595 dot + b;
1596 }
1597 }
1598#endif
1599}
1600
1602 const void *B_prepared,
1603 const float *bias,
1604 float *C,
1605 int M, int N, int K,
1606 int m0, int m1,
1607 int n0, int n1)
1608{
1609 int use_avx512_vnni = 0;
1610#if defined(__AVX512F__) && defined(__AVX512BW__) && \
1611 defined(__AVX512VNNI__)
1612 use_avx512_vnni = 1;
1613#endif
1615 A_q8, B_prepared, bias, C, M, N, K,
1616 m0, m1, n0, n1, use_avx512_vnni);
1617}
1618
1620 const void *A_q8,
1621 const void *B_prepared,
1622 const float *bias,
1623 float *C,
1624 int M, int N, int K)
1625{
1627 A_q8, B_prepared, bias, C, M, N, K,
1628 0, M, 0, N, 1);
1629}
1630
1631void gemm_nt_q6_k_q8_k_prepared(const void *A_q8,
1632 const void *B_prepared,
1633 const float *bias,
1634 float *C,
1635 int M, int N, int K)
1636{
1638 A_q8, B_prepared, bias, C, M, N, K, 0, M, 0, N);
1639}
1640
1641static inline float ck_dot_q6_k_q8_k_fast_or_ref(const block_q6_K *w,
1642 const block_q8_K *x,
1643 int K)
1644{
1646 return dot_q6_k_q8_k_ref(w, x, K);
1647 }
1648#if defined(__AVX2__)
1649 return dot_q6_k_q8_k_avx2(w, x, K);
1650#elif defined(__AVX__)
1651 return dot_q6_k_q8_k_avx(w, x, K);
1652#elif defined(__SSE4_1__)
1653 return dot_q6_k_q8_k_sse(w, x, K);
1654#else
1655 return dot_q6_k_q8_k_ref(w, x, K);
1656#endif
1657}
1658
1659/**
1660 * @brief Compute one C[m0:m1, n0:n1] tile for Q8_K activations x Q6_K weights.
1661 *
1662 * Pure tile math only: no threadpool, no global scheduling, no allocation.
1663 * The orchestrator decides how to split tile jobs across cores.
1664 */
1665void gemm_nt_q6_k_q8_k_tile(const void *A_q8,
1666 const void *B,
1667 const float *bias,
1668 float *C,
1669 int M, int N, int K,
1670 int m0, int m1,
1671 int n0, int n1)
1672{
1673 if (!A_q8 || !B || !C) {
1674 return;
1675 }
1676 if (M <= 0 || N <= 0 || K <= 0 || K % QK_K != 0) {
1677 return;
1678 }
1679 if (m0 < 0) m0 = 0;
1680 if (n0 < 0) n0 = 0;
1681 if (m1 > M) m1 = M;
1682 if (n1 > N) n1 = N;
1683 if (m0 >= m1 || n0 >= n1) {
1684 return;
1685 }
1686
1687 const block_q8_K *A = (const block_q8_K *)A_q8;
1688 const block_q6_K *W = (const block_q6_K *)B;
1689 const int blocks_per_vec = K / QK_K;
1690 const int blocks_per_row = K / QK_K;
1691
1692 for (int n = n0; n < n1; ++n) {
1693 const block_q6_K *w_row = W + (size_t)n * (size_t)blocks_per_row;
1694 const float b = bias ? bias[n] : 0.0f;
1695 for (int m = m0; m < m1; ++m) {
1696 const block_q8_K *a_row = A + (size_t)m * (size_t)blocks_per_vec;
1697 C[(size_t)m * (size_t)N + (size_t)n] = ck_dot_q6_k_q8_k_fast_or_ref(w_row, a_row, K) + b;
1698 }
1699 }
1700}
1701
1702void gemm_nt_q6_k_q8_k_m4_tile(const void *A_q8,
1703 const void *B,
1704 const float *bias,
1705 float *C,
1706 int M, int N, int K,
1707 int m0, int m1,
1708 int n0, int n1)
1709{
1710 if (!A_q8 || !B || !C || M <= 0 || N <= 0 || K <= 0 ||
1711 K % QK_K != 0) {
1712 return;
1713 }
1714 if (m0 < 0) m0 = 0;
1715 if (n0 < 0) n0 = 0;
1716 if (m1 > M) m1 = M;
1717 if (n1 > N) n1 = N;
1718 if (m0 >= m1 || n0 >= n1) return;
1719
1720#if defined(__AVX2__)
1722 const block_q8_K *A = (const block_q8_K *)A_q8;
1723 const block_q6_K *W = (const block_q6_K *)B;
1724 const int blocks_per_vec = K / QK_K;
1725 for (int n = n0; n < n1; ++n) {
1726 const block_q6_K *w_row =
1727 W + (size_t)n * (size_t)blocks_per_vec;
1728 const float b = bias ? bias[n] : 0.0f;
1729 int m = m0;
1730 for (; m + 4 <= m1; m += 4) {
1731 float values[4];
1732 dot_q6_k_q8_k_avx2_m4(
1733 w_row, A + (size_t)m * (size_t)blocks_per_vec,
1734 blocks_per_vec, 4, K, values);
1735 for (int r = 0; r < 4; ++r) {
1736 C[(size_t)(m + r) * (size_t)N + (size_t)n] = values[r] + b;
1737 }
1738 }
1739 if (m < m1) {
1740 float values[4];
1741 const int rows = m1 - m;
1742 dot_q6_k_q8_k_avx2_m4(
1743 w_row, A + (size_t)m * (size_t)blocks_per_vec,
1744 blocks_per_vec, rows, K, values);
1745 for (int r = 0; r < rows; ++r) {
1746 C[(size_t)(m + r) * (size_t)N + (size_t)n] = values[r] + b;
1747 }
1748 }
1749 }
1750 return;
1751 }
1752#endif
1754 A_q8, B, bias, C, M, N, K, m0, m1, n0, n1);
1755}
1756
1757/**
1758 * @brief Experimental single-thread tiled NT GEMM wrapper.
1759 *
1760 * Kept as a separate symbol from gemm_nt_q6_k_q8_k for benchmarks and parity.
1761 * Production prefill should prefer the v8 2D tile scheduler when enabled.
1762 */
1763void gemm_nt_q6_k_q8_k_tiled(const void *A_q8,
1764 const void *B,
1765 const float *bias,
1766 float *C,
1767 int M, int N, int K)
1768{
1769 enum { TILE_M = 8, TILE_N = 16 };
1770 for (int n0 = 0; n0 < N; n0 += TILE_N) {
1771 const int n1 = (n0 + TILE_N < N) ? (n0 + TILE_N) : N;
1772 for (int m0 = 0; m0 < M; m0 += TILE_M) {
1773 const int m1 = (m0 + TILE_M < M) ? (m0 + TILE_M) : M;
1774 gemm_nt_q6_k_q8_k_tile(A_q8, B, bias, C, M, N, K, m0, m1, n0, n1);
1775 }
1776 }
1777}
int ck_strict_parity_enabled(void)
Quantization block structures for weight-only quantization.
#define GGML_FP16_TO_FP32
ck_half ggml_half
#define QK_K
void gemm_nt_q6_k_q8_k_tile(const void *A_q8, const void *B, const float *bias, float *C, int M, int N, int K, int m0, int m1, int n0, int n1)
Compute one C[m0:m1, n0:n1] tile for Q8_K activations x Q6_K weights.
void gemv_q6_k_q8_k_ref(float *y, const void *W, const void *x_q8, int M, int K)
void gemv_q6_k_q8_k_sse(float *y, const void *W, const void *x_q8, int M, int K)
size_t ck_q6_k_prepared_block_size(void)
void gemv_q6_k_q8_k_avx(float *y, const void *W, const void *x_q8, int M, int K)
const char * ck_q6_k_q8_k_provider_name(void)
void gemm_nt_q6_k_q8_k_m4_tile(const void *A_q8, const void *B, const float *bias, float *C, int M, int N, int K, int m0, int m1, int n0, int n1)
static float ck_dot_q6_k_q8_k_fast_or_ref(const block_q6_K *w, const block_q8_K *x, int K)
void gemm_q6_k_q8_k(float *Y, const void *W, const void *X_q8, int M, int N, int K)
GEMM: Y = W @ X^T where W is Q6_K and X is Q8_K.
static void gemm_nt_q6_k_q8_k_prepared_tile_impl(const void *A_q8, const void *B_prepared, const float *bias, float *C, int M, int N, int K, int m0, int m1, int n0, int n1, int use_avx512_vnni)
void gemv_q6_k_q8_k(float *y, const void *W, const void *x_q8, int M, int K)
GEMV: y = W @ x where W is Q6_K and x is Q8_K.
void gemv_q6_k_q8_k_avx2(float *y, const void *W, const void *x_q8, int M, int K)
void vec_dot_q6_k_q8_k(int n, float *s, const void *vx, const void *vy)
Q6_K x Q8_K dot product (single row)
void gemm_nt_q6_k_q8_k_prepared_avx512_vnni(const void *A_q8, const void *B_prepared, const float *bias, float *C, int M, int N, int K)
void gemm_nt_q6_k_q8_k_prepared(const void *A_q8, const void *B_prepared, const float *bias, float *C, int M, int N, int K)
void gemm_nt_q6_k_q8_k_tiled(const void *A_q8, const void *B, const float *bias, float *C, int M, int N, int K)
Experimental single-thread tiled NT GEMM wrapper.
void gemv_q6_k_q8_k_avx512_vbmi(float *y, const void *W, const void *x_q8, int M, int K)
void gemv_q6_k_q8_k_parallel(float *y, const void *W, const void *x_q8, int M, int K, int ith, int nth)
Parallel reference GEMV for Q6_K × Q8_K.
void gemm_nt_q6_k_q8_k_prepared_tile(const void *A_q8, const void *B_prepared, const float *bias, float *C, int M, int N, int K, int m0, int m1, int n0, int n1)
static int ck_q6k_q8k_force_ref(void)
void ck_q6_k_prepare_weight(const void *src, void *dst, int N, int K)
const char * ck_q6_k_prepared_provider_name(void)
void gemm_nt_q6_k_q8_k(const void *A_q8, const void *B, const float *bias, float *C, int M, int N, int K)
NT GEMM: C = A @ B^T where A is Q8_K and B is Q6_K.
void gemv_q6_k_q8_k_avx512(float *y, const void *W, const void *x_q8, int M, int K)
void gemv_q6_k_q8_k_parallel_simd(float *y, const void *W, const void *x_q8, int M, int K, int ith, int nth)
Parallel SIMD GEMV for Q6_K × Q8_K.
static float dot_q6_k_q8_k_ref(const block_q6_K *w, const block_q8_K *x, int K)
Scalar dot product for Q6_K x Q8_K.
#define C(color)
Definition show_config.c:39
uint8_t ql[256/2]
int8_t scales[256/16]
uint8_t qh[256/4]
int8_t qs[256]
int16_t bsums[256/16]