← Back to C-Kernel-Engine Docs Doxygen Source Documentation
 
Loading...
Searching...
No Matches
gemm_kernels_q5_0.c
Go to the documentation of this file.
1/**
2 * @file gemm_kernels_q5_0.c
3 * @brief GEMM/GEMV kernels with Q5_0 quantized weights
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 * Q5_0 Format:
15 * - 32 weights per block
16 * - 1 FP16 scale per block
17 * - Low 4-bits stored like Q4_0 (16 bytes)
18 * - High 1-bit packed separately (4 bytes)
19 * - 22 bytes per 32 weights = 5.5 bits/weight
20 *
21 * Dequantization: w = scale * (q5 - 16)
22 * where q5 = low4bit | (highbit << 4), giving values 0-31, then subtract 16 for signed -16 to +15
23 *
24 * Operations:
25 * Forward: Y = W @ X (W is Q5_0, X and Y are FP32)
26 * Backward: dX = W^T @ dY (gradient w.r.t. input)
27 */
28
29#include <stdint.h>
30#include <stddef.h>
31#include <string.h>
32#include <stdio.h>
33#include "ckernel_quant.h"
34#include "ck_features.h"
35
36/* Include SIMD headers based on available extensions */
37#if defined(__AVX512F__) || defined(__AVX2__) || defined(__AVX__) || defined(__SSE4_1__)
38#include <immintrin.h>
39#endif
40#if defined(__ARM_NEON) || defined(__aarch64__)
41#include <arm_neon.h>
42#endif
43
44/* Forward declarations for dequant functions (defined in dequant_kernels.c) */
45void dequant_q5_0_block(const block_q5_0 *block, float *output);
46void dequant_q5_0_row(const void *src, float *dst, size_t n_elements);
47
48void gemm_nt_q5_0_sse_v2(const float *A,
49 const void *B,
50 const float *bias,
51 float *C,
52 int M, int N, int K);
53
54/* ============================================================================
55 * Forward Pass: GEMV y = W @ x
56 * ============================================================================ */
57
58/**
59 * @brief Matrix-vector multiply with Q5_0 weights (scalar reference)
60 *
61 * @param y Output vector [M]
62 * @param W Weight matrix in Q5_0 format [M x K]
63 * @param x Input vector [K]
64 * @param M Number of output rows
65 * @param K Number of columns (must be multiple of 32)
66 */
67void gemv_q5_0_ref(float *y,
68 const void *W,
69 const float *x,
70 int M, int K)
71{
72 const block_q5_0 *blocks = (const block_q5_0 *)W;
73 const int blocks_per_row = K / QK5_0;
74
75 for (int row = 0; row < M; row++) {
76 float sum = 0.0f;
77
78 for (int b = 0; b < blocks_per_row; b++) {
79 const block_q5_0 *block = &blocks[row * blocks_per_row + b];
80 const float d = CK_FP16_TO_FP32(block->d);
81 const float *xp = &x[b * QK5_0];
82
83 /* Get high bits as 32-bit integer */
84 uint32_t qh;
85 memcpy(&qh, block->qh, sizeof(qh));
86
87 /* llama.cpp Q5_0 layout:
88 * - Weight j uses: low nibble of qs[j], high bit from qh bit j
89 * - Weight j+16 uses: high nibble of qs[j], high bit from qh bit (j+12)
90 * Note: j+12 not j+16 for the high bit of the second weight!
91 */
92 for (int j = 0; j < QK5_0 / 2; j++) {
93 const uint8_t packed = block->qs[j];
94
95 /* Extract nibbles */
96 const int lo = (packed & 0x0F);
97 const int hi = (packed >> 4);
98
99 /* Extract high bits - matches llama.cpp exactly */
100 const int xh_0 = ((qh >> (j + 0)) << 4) & 0x10;
101 const int xh_1 = ((qh >> (j + 12))) & 0x10;
102
103 /* Combine to 5-bit signed value */
104 const int q0 = (lo | xh_0) - 16;
105 const int q1 = (hi | xh_1) - 16;
106
107 /* Weights at indices j and j+16 */
108 sum += d * (float)q0 * xp[j];
109 sum += d * (float)q1 * xp[j + 16];
110 }
111 }
112
113 y[row] = sum;
114 }
115}
116
117#ifdef __AVX512F__
118/**
119 * @brief Matrix-vector multiply with Q5_0 weights (AVX-512)
120 */
121void gemv_q5_0_avx512(float *y,
122 const void *W,
123 const float *x,
124 int M, int K)
125{
126 const block_q5_0 *blocks = (const block_q5_0 *)W;
127 const int blocks_per_row = K / QK5_0;
128 const __m512i offset = _mm512_set1_epi32(16);
129 const __m512i mask_lo = _mm512_set1_epi32(0x0F);
130 const __m512i one = _mm512_set1_epi32(1);
131
132 for (int row = 0; row < M; row++) {
133 __m512 acc = _mm512_setzero_ps();
134
135 for (int b = 0; b < blocks_per_row; b++) {
136 const block_q5_0 *block = &blocks[row * blocks_per_row + b];
137 const __m512 vscale = _mm512_set1_ps(CK_FP16_TO_FP32(block->d));
138 const float *xp = &x[b * QK5_0];
139
140 /* Load high bits */
141 uint32_t qh;
142 memcpy(&qh, block->qh, sizeof(qh));
143
144 /* Load 16 bytes = 32 x 4-bit low weights */
145 __m128i packed = _mm_loadu_si128((const __m128i *)block->qs);
146 __m512i bytes = _mm512_cvtepu8_epi32(packed);
147
148 /* Extract low nibbles */
149 __m512i lo = _mm512_and_epi32(bytes, mask_lo);
150 __m512i hi_shift = _mm512_srli_epi32(bytes, 4);
151
152 /* llama.cpp Q5_0 layout:
153 * - Weights 0-15: high bits from qh bits 0-15
154 * - Weights 16-31: high bits from qh bits 12-27 (j+12 where j=0..15)
155 */
156 /* Build high bit contribution for first 16 weights (indices 0-15) */
157 __m512i qh_lo = _mm512_set_epi32(
158 ((qh >> 15) & 1) << 4, ((qh >> 14) & 1) << 4,
159 ((qh >> 13) & 1) << 4, ((qh >> 12) & 1) << 4,
160 ((qh >> 11) & 1) << 4, ((qh >> 10) & 1) << 4,
161 ((qh >> 9) & 1) << 4, ((qh >> 8) & 1) << 4,
162 ((qh >> 7) & 1) << 4, ((qh >> 6) & 1) << 4,
163 ((qh >> 5) & 1) << 4, ((qh >> 4) & 1) << 4,
164 ((qh >> 3) & 1) << 4, ((qh >> 2) & 1) << 4,
165 ((qh >> 1) & 1) << 4, ((qh >> 0) & 1) << 4
166 );
167
168 /* Build high bit contribution for second 16 weights (indices 16-31)
169 * Scalar code: xh_1 = ((qh >> (j + 12))) & 0x10 extracts bit (j+16)
170 * (since 0x10 = bit 4, position is j+12+4 = j+16)
171 * So weights 16-31 use qh bits 16-31 */
172 __m512i qh_hi = _mm512_set_epi32(
173 ((qh >> 31) & 1) << 4, ((qh >> 30) & 1) << 4,
174 ((qh >> 29) & 1) << 4, ((qh >> 28) & 1) << 4,
175 ((qh >> 27) & 1) << 4, ((qh >> 26) & 1) << 4,
176 ((qh >> 25) & 1) << 4, ((qh >> 24) & 1) << 4,
177 ((qh >> 23) & 1) << 4, ((qh >> 22) & 1) << 4,
178 ((qh >> 21) & 1) << 4, ((qh >> 20) & 1) << 4,
179 ((qh >> 19) & 1) << 4, ((qh >> 18) & 1) << 4,
180 ((qh >> 17) & 1) << 4, ((qh >> 16) & 1) << 4
181 );
182
183 /* Combine low + high bits and subtract offset */
184 __m512i q_lo = _mm512_sub_epi32(_mm512_or_epi32(lo, qh_lo), offset);
185 __m512i q_hi = _mm512_sub_epi32(_mm512_or_epi32(hi_shift, qh_hi), offset);
186
187 /* Dequantize */
188 __m512 w_lo = _mm512_mul_ps(_mm512_cvtepi32_ps(q_lo), vscale);
189 __m512 w_hi = _mm512_mul_ps(_mm512_cvtepi32_ps(q_hi), vscale);
190
191 /* Load sequential input: x[0-15] and x[16-31] */
192 __m512 x_first = _mm512_loadu_ps(&xp[0]); /* x[0..15] */
193 __m512 x_second = _mm512_loadu_ps(&xp[16]); /* x[16..31] */
194
195 acc = _mm512_fmadd_ps(w_lo, x_first, acc);
196 acc = _mm512_fmadd_ps(w_hi, x_second, acc);
197 }
198
199 y[row] = _mm512_reduce_add_ps(acc);
200 }
201}
202#endif
203
204/* ============================================================================
205 * AVX2 Implementation (Haswell+, 256-bit integer operations)
206 *
207 * AVX2 provides true 256-bit integer operations that AVX lacks.
208 * This implementation uses:
209 * - _mm256_cvtepi8_epi32: Sign-extend 8 bytes to 8 int32s (AVX2)
210 * - _mm256_fmadd_ps: Fused multiply-add (FMA3, available with AVX2)
211 * - 256-bit integer shuffles and masks
212 *
213 * Processing: 32 weights per block, 8 at a time with AVX2 registers
214 * ============================================================================ */
215
216#if defined(__AVX2__) && !defined(__AVX512F__)
217
218/* Helper: AVX2 horizontal sum of 8 floats */
219static inline float hsum_avx2(__m256 v) {
220 __m128 lo = _mm256_castps256_ps128(v);
221 __m128 hi = _mm256_extractf128_ps(v, 1);
222 lo = _mm_add_ps(lo, hi); /* 4 floats */
223 __m128 shuf = _mm_shuffle_ps(lo, lo, _MM_SHUFFLE(2, 3, 0, 1));
224 __m128 sums = _mm_add_ps(lo, shuf);
225 shuf = _mm_movehl_ps(shuf, sums);
226 sums = _mm_add_ss(sums, shuf);
227 return _mm_cvtss_f32(sums);
228}
229
230/**
231 * @brief Matrix-vector multiply with Q5_0 weights (AVX2 optimized)
232 *
233 * Uses AVX2 256-bit integer operations for efficient dequantization.
234 * Processes 8 weights at a time with full 256-bit registers.
235 */
236void gemv_q5_0_avx2(float *y,
237 const void *W,
238 const float *x,
239 int M, int K)
240{
241 const block_q5_0 *blocks = (const block_q5_0 *)W;
242 const int blocks_per_row = K / QK5_0; /* QK5_0 = 32 */
243
244 for (int row = 0; row < M; row++) {
245 __m256 acc = _mm256_setzero_ps();
246
247 for (int b = 0; b < blocks_per_row; b++) {
248 const block_q5_0 *block = &blocks[row * blocks_per_row + b];
249 const float d = CK_FP16_TO_FP32(block->d);
250 const __m256 vscale = _mm256_set1_ps(d);
251 const float *xp = &x[b * QK5_0];
252
253 /* Get high bits as 32-bit integer */
254 uint32_t qh;
255 memcpy(&qh, block->qh, sizeof(qh));
256
257 /* Q5_0 layout: 32 weights per block
258 * - Weights 0-15: low nibbles of qs[0-15], high bit from qh[0-15]
259 * - Weights 16-31: high nibbles of qs[0-15], high bit from qh[16-31]
260 *
261 * Process in 4 groups of 8 for AVX2:
262 */
263
264 /* Group 0: weights 0-7 (low nibbles of qs[0-7], high bits qh[0-7]) */
265 {
266 __m128i qs8 = _mm_loadl_epi64((const __m128i *)block->qs);
267 __m128i lo = _mm_and_si128(qs8, _mm_set1_epi8(0x0F));
268
269 /* Build high bits for weights 0-7 */
270 int8_t hb[8];
271 for (int i = 0; i < 8; i++) {
272 hb[i] = ((qh >> i) << 4) & 0x10;
273 }
274 __m128i hi = _mm_loadl_epi64((const __m128i *)hb);
275
276 /* Combine and subtract offset to get signed values */
277 __m128i q5 = _mm_or_si128(lo, hi);
278 __m128i offset = _mm_set1_epi8(16);
279 __m128i q5_signed = _mm_sub_epi8(q5, offset);
280
281 /* Sign-extend to 32-bit and convert to float (AVX2) */
282 __m256i q32 = _mm256_cvtepi8_epi32(q5_signed);
283 __m256 wf = _mm256_cvtepi32_ps(q32);
284 wf = _mm256_mul_ps(wf, vscale);
285
286 /* Load input and accumulate */
287 __m256 xv = _mm256_loadu_ps(&xp[0]);
288 acc = _mm256_fmadd_ps(wf, xv, acc);
289 }
290
291 /* Group 1: weights 8-15 (low nibbles of qs[8-15], high bits qh[8-15]) */
292 {
293 __m128i qs8 = _mm_loadl_epi64((const __m128i *)(block->qs + 8));
294 __m128i lo = _mm_and_si128(qs8, _mm_set1_epi8(0x0F));
295
296 int8_t hb[8];
297 for (int i = 0; i < 8; i++) {
298 hb[i] = ((qh >> (8 + i)) << 4) & 0x10;
299 }
300 __m128i hi = _mm_loadl_epi64((const __m128i *)hb);
301
302 __m128i q5 = _mm_or_si128(lo, hi);
303 __m128i offset = _mm_set1_epi8(16);
304 __m128i q5_signed = _mm_sub_epi8(q5, offset);
305
306 __m256i q32 = _mm256_cvtepi8_epi32(q5_signed);
307 __m256 wf = _mm256_cvtepi32_ps(q32);
308 wf = _mm256_mul_ps(wf, vscale);
309
310 __m256 xv = _mm256_loadu_ps(&xp[8]);
311 acc = _mm256_fmadd_ps(wf, xv, acc);
312 }
313
314 /* Group 2: weights 16-23 (high nibbles of qs[0-7], high bits qh[16-23]) */
315 {
316 __m128i qs8 = _mm_loadl_epi64((const __m128i *)block->qs);
317 __m128i hi_nib = _mm_and_si128(_mm_srli_epi16(qs8, 4), _mm_set1_epi8(0x0F));
318
319 /* High bits for weights 16-23 come from qh bits 16-23 */
320 int8_t hb[8];
321 for (int i = 0; i < 8; i++) {
322 hb[i] = ((qh >> (16 + i)) & 1) << 4;
323 }
324 __m128i hi = _mm_loadl_epi64((const __m128i *)hb);
325
326 __m128i q5 = _mm_or_si128(hi_nib, hi);
327 __m128i offset = _mm_set1_epi8(16);
328 __m128i q5_signed = _mm_sub_epi8(q5, offset);
329
330 __m256i q32 = _mm256_cvtepi8_epi32(q5_signed);
331 __m256 wf = _mm256_cvtepi32_ps(q32);
332 wf = _mm256_mul_ps(wf, vscale);
333
334 __m256 xv = _mm256_loadu_ps(&xp[16]);
335 acc = _mm256_fmadd_ps(wf, xv, acc);
336 }
337
338 /* Group 3: weights 24-31 (high nibbles of qs[8-15], high bits qh[24-31]) */
339 {
340 __m128i qs8 = _mm_loadl_epi64((const __m128i *)(block->qs + 8));
341 __m128i hi_nib = _mm_and_si128(_mm_srli_epi16(qs8, 4), _mm_set1_epi8(0x0F));
342
343 int8_t hb[8];
344 for (int i = 0; i < 8; i++) {
345 hb[i] = ((qh >> (24 + i)) & 1) << 4;
346 }
347 __m128i hi = _mm_loadl_epi64((const __m128i *)hb);
348
349 __m128i q5 = _mm_or_si128(hi_nib, hi);
350 __m128i offset = _mm_set1_epi8(16);
351 __m128i q5_signed = _mm_sub_epi8(q5, offset);
352
353 __m256i q32 = _mm256_cvtepi8_epi32(q5_signed);
354 __m256 wf = _mm256_cvtepi32_ps(q32);
355 wf = _mm256_mul_ps(wf, vscale);
356
357 __m256 xv = _mm256_loadu_ps(&xp[24]);
358 acc = _mm256_fmadd_ps(wf, xv, acc);
359 }
360 }
361
362 y[row] = hsum_avx2(acc);
363 }
364}
365#endif /* __AVX2__ && !__AVX512F__ */
366
367/* ============================================================================
368 * AVX Implementation with True SIMD Dequantization
369 *
370 * Q5_0 format: 32 weights per block
371 * - d: FP16 scale
372 * - qh: 4 bytes (32 high bits, one per weight)
373 * - qs: 16 bytes (low 4 bits, packed as pairs)
374 * - Dequant: w = d * ((lo | (highbit << 4)) - 16)
375 *
376 * This uses SSE for integer unpacking (Ivy Bridge doesn't have AVX2 for
377 * 256-bit integer ops) and AVX for float accumulation.
378 *
379 * Key optimization: Instead of scalar dequant, we use SIMD to:
380 * 1. Extract nibbles to bytes using SSE shuffle/shift
381 * 2. Combine with high bits using SSE or/and
382 * 3. Convert to float and scale
383 * ============================================================================ */
384
385#if defined(__AVX__) && !defined(__AVX2__) && !defined(__AVX512F__)
386
387/* Helper: Extract low nibbles from 16 packed bytes to 16 bytes */
388static inline __m128i extract_low_nibbles(__m128i packed) {
389 return _mm_and_si128(packed, _mm_set1_epi8(0x0F));
390}
391
392/* Helper: Extract high nibbles from 16 packed bytes to 16 bytes */
393static inline __m128i extract_high_nibbles(__m128i packed) {
394 return _mm_and_si128(_mm_srli_epi16(packed, 4), _mm_set1_epi8(0x0F));
395}
396
397/* Helper: SSE horizontal sum of 4 floats */
398static inline float hsum_sse(__m128 v) {
399 __m128 shuf = _mm_shuffle_ps(v, v, _MM_SHUFFLE(2, 3, 0, 1));
400 __m128 sums = _mm_add_ps(v, shuf);
401 shuf = _mm_movehl_ps(shuf, sums);
402 sums = _mm_add_ss(sums, shuf);
403 return _mm_cvtss_f32(sums);
404}
405
406/* Helper: SSE dot product of 8 int8 values with 8 float values */
407static inline float dot_int8_float8_sse(__m128i q8_lo, const float *x, float scale) {
408 /* Sign-extend 8 int8s to 8 int32s (in two steps) */
409 __m128i lo16 = _mm_cvtepi8_epi16(q8_lo); /* 8 int8 -> 8 int16 */
410 __m128i lo32_0 = _mm_cvtepi16_epi32(lo16); /* low 4 int16 -> 4 int32 */
411 __m128i lo32_1 = _mm_cvtepi16_epi32(_mm_srli_si128(lo16, 8)); /* high 4 int16 -> 4 int32 */
412
413 /* Convert to float */
414 __m128 w0 = _mm_cvtepi32_ps(lo32_0);
415 __m128 w1 = _mm_cvtepi32_ps(lo32_1);
416
417 /* Scale */
418 __m128 vscale = _mm_set1_ps(scale);
419 w0 = _mm_mul_ps(w0, vscale);
420 w1 = _mm_mul_ps(w1, vscale);
421
422 /* Load input and multiply */
423 __m128 x0 = _mm_loadu_ps(x);
424 __m128 x1 = _mm_loadu_ps(x + 4);
425
426 __m128 prod0 = _mm_mul_ps(w0, x0);
427 __m128 prod1 = _mm_mul_ps(w1, x1);
428
429 /* Sum */
430 __m128 sum = _mm_add_ps(prod0, prod1);
431 return hsum_sse(sum);
432}
433
434/**
435 * @brief Matrix-vector multiply with Q5_0 weights (AVX + SSE optimized)
436 *
437 * Uses SSE for integer dequantization, AVX for float accumulation.
438 * ~3-5x faster than scalar reference on Ivy Bridge.
439 */
440void gemv_q5_0_avx(float *y,
441 const void *W,
442 const float *x,
443 int M, int K)
444{
445 const block_q5_0 *blocks = (const block_q5_0 *)W;
446 const int blocks_per_row = K / QK5_0; /* QK5_0 = 32 */
447
448 const __m128i mask_0f = _mm_set1_epi8(0x0F);
449 const __m128i mask_10 = _mm_set1_epi8(0x10);
450
451 for (int row = 0; row < M; row++) {
452 float sum = 0.0f;
453
454 for (int b = 0; b < blocks_per_row; b++) {
455 const block_q5_0 *block = &blocks[row * blocks_per_row + b];
456 const float d = CK_FP16_TO_FP32(block->d);
457 const float *xp = &x[b * QK5_0];
458
459 /* Load 16 packed bytes (32 nibbles) */
460 __m128i qs = _mm_loadu_si128((const __m128i *)block->qs);
461
462 /* Extract low and high nibbles */
463 __m128i lo_nibbles = _mm_and_si128(qs, mask_0f); /* 16 low nibbles */
464 __m128i hi_nibbles = _mm_and_si128(_mm_srli_epi16(qs, 4), mask_0f); /* 16 high nibbles */
465
466 /* Get high bits as 32-bit integer */
467 uint32_t qh;
468 memcpy(&qh, block->qh, sizeof(qh));
469
470 /* Q5_0 layout: weight j uses qs[j/2] nibble (low if j<16, high if j>=16)
471 * plus high bit from qh:
472 * - weights 0-15: low nibbles of qs[0-15], high bit at qh[0-15]
473 * - weights 16-31: high nibbles of qs[0-15], high bit at qh[12-27]
474 *
475 * For efficiency, we process 32 weights in 4 groups of 8:
476 */
477
478 /* Group 0: weights 0-7 (low nibbles of qs[0-7], high bits from qh[0-7]) */
479 {
480 uint8_t w8[8];
481 for (int i = 0; i < 8; i++) {
482 int lo = block->qs[i] & 0x0F;
483 int hb = ((qh >> i) << 4) & 0x10;
484 w8[i] = (lo | hb) - 16; /* Signed -16 to +15 */
485 }
486 __m128i q8 = _mm_loadl_epi64((const __m128i *)w8);
487 sum += dot_int8_float8_sse(q8, &xp[0], d);
488 }
489
490 /* Group 1: weights 8-15 (low nibbles of qs[8-15], high bits from qh[8-15]) */
491 {
492 uint8_t w8[8];
493 for (int i = 0; i < 8; i++) {
494 int lo = block->qs[8 + i] & 0x0F;
495 int hb = ((qh >> (8 + i)) << 4) & 0x10;
496 w8[i] = (lo | hb) - 16;
497 }
498 __m128i q8 = _mm_loadl_epi64((const __m128i *)w8);
499 sum += dot_int8_float8_sse(q8, &xp[8], d);
500 }
501
502 /* Group 2: weights 16-23 (high nibbles of qs[0-7], high bits from qh[12-19]) */
503 {
504 uint8_t w8[8];
505 for (int i = 0; i < 8; i++) {
506 int hi = block->qs[i] >> 4;
507 int hb = (qh >> (12 + i)) & 0x10;
508 w8[i] = (hi | hb) - 16;
509 }
510 __m128i q8 = _mm_loadl_epi64((const __m128i *)w8);
511 sum += dot_int8_float8_sse(q8, &xp[16], d);
512 }
513
514 /* Group 3: weights 24-31 (high nibbles of qs[8-15], high bits from qh[20-27]) */
515 {
516 uint8_t w8[8];
517 for (int i = 0; i < 8; i++) {
518 int hi = block->qs[8 + i] >> 4;
519 int hb = (qh >> (20 + i)) & 0x10;
520 w8[i] = (hi | hb) - 16;
521 }
522 __m128i q8 = _mm_loadl_epi64((const __m128i *)w8);
523 sum += dot_int8_float8_sse(q8, &xp[24], d);
524 }
525 }
526
527 y[row] = sum;
528 }
529}
530#endif /* __AVX__ && !__AVX512F__ */
531
532/**
533 * @brief Auto-dispatch GEMV for Q5_0 weights based on CPU features
534 *
535 * Dispatch priority (best available):
536 * 1. AVX-512 (512-bit vectors) - Intel Skylake-X+
537 * 2. AVX2+FMA (256-bit vectors) - Intel Haswell+
538 * 3. AVX (256-bit vectors) - Intel Sandy Bridge+
539 * 4. SSE4.1 (128-bit vectors) - Intel Nehalem+
540 * 5. Reference (scalar) - Fallback
541 *
542 * Uses ck_features.h for standardized feature detection.
543 *
544 * @param y Output vector [M]
545 * @param W Weight matrix in Q5_0 format [M x K]
546 * @param x Input vector [K]
547 * @param M Number of output rows
548 * @param K Number of input columns (hidden dimension)
549 */
550void gemv_q5_0(float *y,
551 const void *W,
552 const float *x,
553 int M, int K)
554{
555// Dispatch order: AVX512 > AVX2 > AVX > SSE > ref
556#if defined(__AVX512F__)
557 gemv_q5_0_avx512(y, W, x, M, K);
558#elif defined(__AVX2__)
559 gemv_q5_0_avx2(y, W, x, M, K);
560#elif defined(__AVX__)
561 gemv_q5_0_avx(y, W, x, M, K);
562#elif defined(__SSE4_1__)
563 gemv_q5_0_ref(y, W, x, M, K);
564#else
565 gemv_q5_0_ref(y, W, x, M, K);
566#endif
567}
568
569/* ============================================================================
570 * PARALLEL VERSIONS (for parallel orchestration)
571 *
572 * These receive ith (thread index) and nth (total threads) from orchestration.
573 * OpenMP lives in orchestration layer, NOT here.
574 * ============================================================================ */
575
576/**
577 * @brief Parallel reference GEMV for Q5_0 × FP32
578 */
579void gemv_q5_0_parallel(float *y,
580 const void *W,
581 const float *x,
582 int M, int K,
583 int ith, int nth)
584{
585 if (!y || !W || !x || M <= 0 || K <= 0) return;
586 if (ith < 0 || nth <= 0 || ith >= nth) return;
587
588 const int dr = (M + nth - 1) / nth;
589 const int r0 = dr * ith;
590 const int r1 = (r0 + dr < M) ? (r0 + dr) : M;
591
592 if (r0 >= M) return;
593
594 const block_q5_0 *blocks = (const block_q5_0 *)W;
595 const int blocks_per_row = K / QK5_0;
596
597 for (int row = r0; row < r1; row++) {
598 float sum = 0.0f;
599 for (int b = 0; b < blocks_per_row; b++) {
600 const block_q5_0 *block = &blocks[row * blocks_per_row + b];
601 const float d = CK_FP16_TO_FP32(block->d);
602 const float *xp = &x[b * QK5_0];
603
604 uint32_t qh;
605 memcpy(&qh, block->qh, sizeof(qh));
606
607 for (int j = 0; j < QK5_0 / 2; j++) {
608 const uint8_t packed = block->qs[j];
609 const int lo = (packed & 0x0F);
610 const int hi = (packed >> 4);
611 const int xh_0 = ((qh >> (j + 0)) << 4) & 0x10;
612 const int xh_1 = ((qh >> (j + 12))) & 0x10;
613 const int w0 = (lo | xh_0) - 16;
614 const int w1 = (hi | xh_1) - 16;
615 sum += d * (w0 * xp[j] + w1 * xp[j + QK5_0/2]);
616 }
617 }
618 y[row] = sum;
619 }
620}
621
622/**
623 * @brief Parallel SIMD GEMV for Q5_0 × FP32 with prefetching
624 */
625void gemv_q5_0_parallel_simd(float *y,
626 const void *W,
627 const float *x,
628 int M, int K,
629 int ith, int nth)
630{
631 if (!y || !W || !x || M <= 0 || K <= 0) return;
632 if (ith < 0 || nth <= 0 || ith >= nth) return;
633
634 const int dr = (M + nth - 1) / nth;
635 const int r0 = dr * ith;
636 const int r1 = (r0 + dr < M) ? (r0 + dr) : M;
637
638 if (r0 >= M) return;
639
640#if defined(__AVX__) || defined(__SSE4_1__)
641 const block_q5_0 *blocks = (const block_q5_0 *)W;
642 const int blocks_per_row = K / QK5_0;
643 /* Prefetch first few rows */
644 const int PREFETCH_ROWS = 4;
645 for (int p = 0; p < PREFETCH_ROWS && r0 + p < r1; ++p) {
646 const char *row_ptr = (const char *)(blocks + (r0 + p) * blocks_per_row);
647 _mm_prefetch(row_ptr, _MM_HINT_T0);
648 _mm_prefetch(row_ptr + 64, _MM_HINT_T0);
649 }
650
651 for (int row = r0; row < r1; ++row) {
652 /* Prefetch rows ahead */
653 if (row + PREFETCH_ROWS < r1) {
654 const char *prefetch_ptr = (const char *)(blocks + (row + PREFETCH_ROWS) * blocks_per_row);
655 _mm_prefetch(prefetch_ptr, _MM_HINT_T0);
656 _mm_prefetch(prefetch_ptr + 64, _MM_HINT_T0);
657 }
658
659 /* Use SIMD dot product for this row */
660#if defined(__AVX512F__)
661 /* Call single-row AVX512 implementation */
662 gemv_q5_0_avx512(&y[row], (const char *)blocks + row * blocks_per_row * sizeof(block_q5_0), x, 1, K);
663#elif defined(__AVX2__)
664 gemv_q5_0_avx2(&y[row], (const char *)blocks + row * blocks_per_row * sizeof(block_q5_0), x, 1, K);
665#elif defined(__AVX__)
666 gemv_q5_0_avx(&y[row], (const char *)blocks + row * blocks_per_row * sizeof(block_q5_0), x, 1, K);
667#else
668 gemv_q5_0_ref(&y[row], (const char *)blocks + row * blocks_per_row * sizeof(block_q5_0), x, 1, K);
669#endif
670 }
671#else
672 /* Fallback to reference parallel */
673 gemv_q5_0_parallel(y, W, x, M, K, ith, nth);
674#endif
675}
676
677/* ============================================================================
678 * Forward Pass: GEMM Y = W @ X
679 * ============================================================================ */
680
681/**
682 * @brief Matrix-matrix multiply with Q5_0 weights
683 */
684void gemm_q5_0(float *Y,
685 const void *W,
686 const float *X,
687 int M, int N, int K)
688{
689 for (int n = 0; n < N; n++) {
690 gemv_q5_0(&Y[n * M], W, &X[n * K], M, K);
691 }
692}
693
694/* ============================================================================
695 * Backward Pass: Gradient w.r.t. Input
696 * ============================================================================ */
697
698/**
699 * @brief Backward pass: compute input gradient
700 *
701 * @param dX Output gradient w.r.t. input [K]
702 * @param W Weight matrix in Q5_0 format [M x K]
703 * @param dY Gradient w.r.t. output [M]
704 * @param M Number of output rows
705 * @param K Number of columns (input dimension)
706 */
707void gemv_q5_0_backward_ref(float *dX,
708 const void *W,
709 const float *dY,
710 int M, int K)
711{
712 const block_q5_0 *blocks = (const block_q5_0 *)W;
713 const int blocks_per_row = K / QK5_0;
714
715 /* Zero output gradient */
716 memset(dX, 0, K * sizeof(float));
717
718 /* Accumulate: dX += W^T @ dY */
719 for (int row = 0; row < M; row++) {
720 const float dy = dY[row];
721
722 for (int b = 0; b < blocks_per_row; b++) {
723 const block_q5_0 *block = &blocks[row * blocks_per_row + b];
724 const float d = CK_FP16_TO_FP32(block->d);
725 float *dxp = &dX[b * QK5_0];
726
727 /* Get high bits */
728 uint32_t qh;
729 memcpy(&qh, block->qh, sizeof(qh));
730
731 /* llama.cpp Q5_0 layout - note j+12 for second weight high bit */
732 for (int j = 0; j < QK5_0 / 2; j++) {
733 const uint8_t packed = block->qs[j];
734
735 /* Extract and reconstruct 5-bit values */
736 const int lo = (packed & 0x0F);
737 const int hi = (packed >> 4);
738 const int xh_0 = ((qh >> (j + 0)) << 4) & 0x10;
739 const int xh_1 = ((qh >> (j + 12))) & 0x10;
740 const int q0 = (lo | xh_0) - 16;
741 const int q1 = (hi | xh_1) - 16;
742
743 dxp[j] += d * (float)q0 * dy;
744 dxp[j + 16] += d * (float)q1 * dy;
745 }
746 }
747 }
748}
749
750/**
751 * @brief Auto-dispatch backward
752 */
753void gemv_q5_0_backward(float *dX,
754 const void *W,
755 const float *dY,
756 int M, int K)
757{
758 gemv_q5_0_backward_ref(dX, W, dY, M, K);
759}
760
761/**
762 * @brief Batched backward pass
763 */
764void gemm_q5_0_backward(float *dX,
765 const void *W,
766 const float *dY,
767 int M, int N, int K)
768{
769 for (int n = 0; n < N; n++) {
770 gemv_q5_0_backward(&dX[n * K], W, &dY[n * M], M, K);
771 }
772}
773
774/* ============================================================================
775 * GEMM NT (Non-Transpose A, Transpose B) - C = A @ B^T
776 * For inference: A is activations [M x K], B is weights [N x K]
777 * ============================================================================ */
778
779/**
780 * @brief GEMM with transposed Q5_0 weights: C = A @ B^T
781 *
782 * @param A Input activations [M x K], row-major FP32
783 * @param B Weight matrix in Q5_0 format [N x K], row-major quantized
784 * @param bias Optional bias [N], NULL if not used
785 * @param C Output [M x N], row-major FP32
786 * @param M Batch size (number of tokens)
787 * @param N Output dimension (number of rows in B)
788 * @param K Input dimension
789 */
790void gemm_nt_q5_0_ref(const float *A,
791 const void *B,
792 const float *bias,
793 float *C,
794 int M, int N, int K)
795{
796 const block_q5_0 *blocks = (const block_q5_0 *)B;
797 const int blocks_per_row = K / QK5_0;
798
799 for (int m = 0; m < M; m++) {
800 const float *a_row = &A[m * K];
801
802 for (int n = 0; n < N; n++) {
803 float sum = 0.0f;
804
805 for (int b = 0; b < blocks_per_row; b++) {
806 const block_q5_0 *block = &blocks[n * blocks_per_row + b];
807 const float d = CK_FP16_TO_FP32(block->d);
808 const float *ap = &a_row[b * QK5_0];
809
810 uint32_t qh;
811 memcpy(&qh, block->qh, sizeof(qh));
812
813 /* llama.cpp Q5_0 layout - note j+12 for second weight high bit */
814 for (int j = 0; j < QK5_0 / 2; j++) {
815 const uint8_t packed = block->qs[j];
816 const int lo = (packed & 0x0F);
817 const int hi = (packed >> 4);
818 const int xh_0 = ((qh >> (j + 0)) << 4) & 0x10;
819 const int xh_1 = ((qh >> (j + 12))) & 0x10;
820 const int q0 = (lo | xh_0) - 16;
821 const int q1 = (hi | xh_1) - 16;
822
823 sum += d * (float)q0 * ap[j];
824 sum += d * (float)q1 * ap[j + 16];
825 }
826 }
827
828 C[m * N + n] = sum + (bias ? bias[n] : 0.0f);
829 }
830 }
831}
832
833void gemm_nt_q5_0(const float *A,
834 const void *B,
835 const float *bias,
836 float *C,
837 int M, int N, int K)
838{
839 /* For decode (M=1), use direct GEMV which has AVX optimization */
840 if (M == 1) {
841 /* gemm_q5_0 expects column-major output, but we need row-major
842 * So we call gemv_q5_0 directly for each output element */
843 gemv_q5_0(C, B, A, N, K);
844 if (bias) {
845 for (int n = 0; n < N; n++) {
846 C[n] += bias[n];
847 }
848 }
849 return;
850 }
851
852 /* For prefill (M>1), use GEMM which dispatches to GEMV with AVX/AVX512 */
853 /* gemm_q5_0 produces Y as [batch x M_out]. Here:
854 * batch = M (tokens)
855 * M_out = N (output channels) */
856 gemm_q5_0(C, B, A, /*M_out=*/N, /*N_batch=*/M, K);
857
858 if (bias) {
859 for (int m = 0; m < M; m++) {
860 float *row = C + (size_t)m * (size_t)N;
861 for (int n = 0; n < N; n++) {
862 row[n] += bias[n];
863 }
864 }
865 }
866}
867
868/* ============================================================================
869 * Dot Product Utility
870 * ============================================================================ */
871
872float dot_q5_0(const void *w_q5_0, const float *x, int K)
873{
874 float result;
875 gemv_q5_0(&result, w_q5_0, x, 1, K);
876 return result;
877}
878
879/* ============================================================================
880 * Quantized Dot Product: Q5_0 x Q8_0
881 *
882 * This matches llama.cpp's ggml_vec_dot_q5_0_q8_0 exactly.
883 * Input is pre-quantized to Q8_0 format, enabling integer dot products.
884 * Result: sum_blocks( (d_w * d_x) * sum_weights( w5 * x8 ) )
885 *
886 * Key difference from gemv_q5_0:
887 * - gemv_q5_0: Takes FP32 input, dequantizes weights to FP32, FP32 dot
888 * - vec_dot_q5_0_q8_0: Takes Q8_0 input, does integer dot, scales at end
889 *
890 * The quantized path is faster and matches llama.cpp for parity testing.
891 * ============================================================================ */
892
893/**
894 * @brief Quantized dot product: Q5_0 weights x Q8_0 input (scalar reference)
895 *
896 * @param n Number of elements (must be multiple of 32)
897 * @param s Output: scalar dot product result
898 * @param vx Q5_0 quantized weights
899 * @param vy Q8_0 quantized input
900 */
901void vec_dot_q5_0_q8_0_ref(int n, float *s, const void *vx, const void *vy)
902{
903 const int qk = QK5_0; /* 32 */
904 const int nb = n / qk;
905
906 const block_q5_0 *x = (const block_q5_0 *)vx;
907 const block_q8_0 *y = (const block_q8_0 *)vy;
908
909 float sumf = 0.0f;
910
911 for (int ib = 0; ib < nb; ib++) {
912 /* Load high bits for this block */
913 uint32_t qh;
914 memcpy(&qh, x[ib].qh, sizeof(qh));
915
916 int sumi0 = 0;
917 int sumi1 = 0;
918
919 for (int j = 0; j < qk / 2; j++) {
920 /* Extract high bits - matches llama.cpp exactly */
921 const uint8_t xh_0 = ((qh & (1u << (j + 0))) >> (j + 0)) << 4;
922 const uint8_t xh_1 = ((qh & (1u << (j + 16))) >> (j + 12));
923
924 /* Reconstruct 5-bit signed values (-16 to +15) */
925 const int32_t x0 = (int8_t)(((x[ib].qs[j] & 0x0F) | xh_0) - 16);
926 const int32_t x1 = (int8_t)(((x[ib].qs[j] >> 4) | xh_1) - 16);
927
928 /* Integer dot product with Q8_0 values */
929 sumi0 += x0 * y[ib].qs[j];
930 sumi1 += x1 * y[ib].qs[j + qk / 2];
931 }
932
933 int sumi = sumi0 + sumi1;
934 sumf += (CK_FP16_TO_FP32(x[ib].d) * CK_FP16_TO_FP32(y[ib].d)) * sumi;
935 }
936
937 *s = sumf;
938}
939
940#if defined(__ARM_NEON) || defined(__aarch64__)
941static inline int32_t ck_hsum_s32x4(int32x4_t v)
942{
943 int32_t lanes[4];
944 vst1q_s32(lanes, v);
945 return lanes[0] + lanes[1] + lanes[2] + lanes[3];
946}
947
948void vec_dot_q5_0_q8_0_neon(int n, float *s, const void *vx, const void *vy)
949{
950 const int qk = QK5_0;
951 const int nb = n / qk;
952
953 const block_q5_0 *x = (const block_q5_0 *)vx;
954 const block_q8_0 *y = (const block_q8_0 *)vy;
955
956 float sumf = 0.0f;
957
958 for (int ib = 0; ib < nb; ib++) {
959 uint32_t qh;
960 memcpy(&qh, x[ib].qh, sizeof(qh));
961
962 int8_t wvals[QK5_0];
963 for (int j = 0; j < qk / 2; j++) {
964 const uint8_t packed = x[ib].qs[j];
965 const uint8_t xh_0 = ((qh >> (j + 0)) & 1u) << 4;
966 const uint8_t xh_1 = ((qh >> (j + 16)) & 1u) << 4;
967
968 wvals[j] = (int8_t)(((packed & 0x0F) | xh_0) - 16);
969 wvals[j + qk / 2] = (int8_t)(((packed >> 4) | xh_1) - 16);
970 }
971
972 const int8x16_t w0 = vld1q_s8(&wvals[0]);
973 const int8x16_t w1 = vld1q_s8(&wvals[16]);
974 const int8x16_t x0 = vld1q_s8(&y[ib].qs[0]);
975 const int8x16_t x1 = vld1q_s8(&y[ib].qs[16]);
976
977 int32x4_t acc = vdupq_n_s32(0);
978
979 int16x8_t p0 = vmull_s8(vget_low_s8(w0), vget_low_s8(x0));
980 int16x8_t p1 = vmull_s8(vget_high_s8(w0), vget_high_s8(x0));
981 int16x8_t p2 = vmull_s8(vget_low_s8(w1), vget_low_s8(x1));
982 int16x8_t p3 = vmull_s8(vget_high_s8(w1), vget_high_s8(x1));
983
984 acc = vaddq_s32(acc, vpaddlq_s16(p0));
985 acc = vaddq_s32(acc, vpaddlq_s16(p1));
986 acc = vaddq_s32(acc, vpaddlq_s16(p2));
987 acc = vaddq_s32(acc, vpaddlq_s16(p3));
988
989 const float d = CK_FP16_TO_FP32(x[ib].d) * CK_FP16_TO_FP32(y[ib].d);
990 sumf += d * (float)ck_hsum_s32x4(acc);
991 }
992
993 *s = sumf;
994}
995#endif
996
997#ifdef __AVX512F__
998/**
999 * @brief Quantized dot product Q5_0 x Q8_0 (AVX-512)
1000 */
1001void vec_dot_q5_0_q8_0_avx512(int n, float *s, const void *vx, const void *vy)
1002{
1003 const int qk = QK5_0;
1004 const int nb = n / qk;
1005
1006 const block_q5_0 *x = (const block_q5_0 *)vx;
1007 const block_q8_0 *y = (const block_q8_0 *)vy;
1008
1009 float sumf = 0.0f;
1010
1011 for (int ib = 0; ib < nb; ib++) {
1012 const float d = CK_FP16_TO_FP32(x[ib].d) * CK_FP16_TO_FP32(y[ib].d);
1013
1014 /* Load high bits */
1015 uint32_t qh;
1016 memcpy(&qh, x[ib].qh, sizeof(qh));
1017
1018 /* Load 16 packed bytes (32 nibbles) */
1019 __m128i qs = _mm_loadu_si128((const __m128i *)x[ib].qs);
1020
1021 /* Process first 16 weights (low nibbles, high bits 0-15) */
1022 __m512i lo_nibbles = _mm512_cvtepu8_epi32(qs);
1023 lo_nibbles = _mm512_and_epi32(lo_nibbles, _mm512_set1_epi32(0x0F));
1024
1025 /* Build high bit contribution for first 16 weights */
1026 __m512i qh_lo = _mm512_set_epi32(
1027 ((qh >> 15) & 1) << 4, ((qh >> 14) & 1) << 4,
1028 ((qh >> 13) & 1) << 4, ((qh >> 12) & 1) << 4,
1029 ((qh >> 11) & 1) << 4, ((qh >> 10) & 1) << 4,
1030 ((qh >> 9) & 1) << 4, ((qh >> 8) & 1) << 4,
1031 ((qh >> 7) & 1) << 4, ((qh >> 6) & 1) << 4,
1032 ((qh >> 5) & 1) << 4, ((qh >> 4) & 1) << 4,
1033 ((qh >> 3) & 1) << 4, ((qh >> 2) & 1) << 4,
1034 ((qh >> 1) & 1) << 4, ((qh >> 0) & 1) << 4
1035 );
1036
1037 /* Combine and subtract 16 to get signed values */
1038 __m512i q5_lo = _mm512_sub_epi32(_mm512_or_epi32(lo_nibbles, qh_lo),
1039 _mm512_set1_epi32(16));
1040
1041 /* Load Q8_0 values for first 16 */
1042 __m128i y8_lo = _mm_loadu_si128((const __m128i *)&y[ib].qs[0]);
1043 __m512i y32_lo = _mm512_cvtepi8_epi32(y8_lo);
1044
1045 /* Integer multiply and accumulate */
1046 __m512i prod_lo = _mm512_mullo_epi32(q5_lo, y32_lo);
1047
1048 /* Process second 16 weights (high nibbles, high bits 16-31 via j+12 mapping) */
1049 __m512i hi_nibbles = _mm512_cvtepu8_epi32(qs);
1050 hi_nibbles = _mm512_srli_epi32(hi_nibbles, 4);
1051
1052 /* Build high bit contribution for second 16 weights (bits 16-31) */
1053 __m512i qh_hi = _mm512_set_epi32(
1054 ((qh >> 31) & 1) << 4, ((qh >> 30) & 1) << 4,
1055 ((qh >> 29) & 1) << 4, ((qh >> 28) & 1) << 4,
1056 ((qh >> 27) & 1) << 4, ((qh >> 26) & 1) << 4,
1057 ((qh >> 25) & 1) << 4, ((qh >> 24) & 1) << 4,
1058 ((qh >> 23) & 1) << 4, ((qh >> 22) & 1) << 4,
1059 ((qh >> 21) & 1) << 4, ((qh >> 20) & 1) << 4,
1060 ((qh >> 19) & 1) << 4, ((qh >> 18) & 1) << 4,
1061 ((qh >> 17) & 1) << 4, ((qh >> 16) & 1) << 4
1062 );
1063
1064 __m512i q5_hi = _mm512_sub_epi32(_mm512_or_epi32(hi_nibbles, qh_hi),
1065 _mm512_set1_epi32(16));
1066
1067 /* Load Q8_0 values for second 16 */
1068 __m128i y8_hi = _mm_loadu_si128((const __m128i *)&y[ib].qs[16]);
1069 __m512i y32_hi = _mm512_cvtepi8_epi32(y8_hi);
1070
1071 __m512i prod_hi = _mm512_mullo_epi32(q5_hi, y32_hi);
1072
1073 /* Sum all products */
1074 int sumi = _mm512_reduce_add_epi32(_mm512_add_epi32(prod_lo, prod_hi));
1075
1076 /* Scale and accumulate - use scalar to avoid 16x broadcast bug */
1077 sumf += d * (float)sumi;
1078 }
1079
1080 *s = sumf;
1081}
1082#endif
1083
1084#if defined(__AVX2__)
1085static inline __m256i bytes_from_bits_32_avx(const uint8_t *qh);
1086static inline __m256i bytes_from_nibbles_32_avx(const uint8_t *qs);
1087static inline __m256 mul_sum_i8_pairs_float_avx(const __m256i x, const __m256i y);
1088static inline float hsum_float_8_avx(const __m256 x);
1089
1090/**
1091 * @brief Quantized dot product Q5_0 x Q8_0 (AVX2/FMA)
1092 *
1093 * Mirrors the local llama.cpp AVX2 structure on machines like this one: unpack
1094 * Q5_0 into signed int8 lanes, multiply against Q8_0 int8 lanes, and fuse the
1095 * per-block scale into the FP32 accumulator.
1096 */
1097void vec_dot_q5_0_q8_0_avx2(int n, float *s, const void *vx, const void *vy)
1098{
1099 const int qk = QK5_0;
1100 const int nb = n / qk;
1101
1102 const block_q5_0 *x = (const block_q5_0 *)vx;
1103 const block_q8_0 *y = (const block_q8_0 *)vy;
1104
1105 __m256 acc = _mm256_setzero_ps();
1106
1107 for (int ib = 0; ib < nb; ++ib) {
1108 const __m256 d =
1109 _mm256_set1_ps(CK_FP16_TO_FP32(x[ib].d) * CK_FP16_TO_FP32(y[ib].d));
1110
1111 __m256i qx = bytes_from_nibbles_32_avx(x[ib].qs);
1112 __m256i bxhi = bytes_from_bits_32_avx(x[ib].qh);
1113 bxhi = _mm256_andnot_si256(bxhi, _mm256_set1_epi8((char)0xF0));
1114 qx = _mm256_or_si256(qx, bxhi);
1115
1116 const __m256i qy = _mm256_loadu_si256((const __m256i *)y[ib].qs);
1117 const __m256 q = mul_sum_i8_pairs_float_avx(qx, qy);
1118#if defined(__FMA__)
1119 acc = _mm256_fmadd_ps(d, q, acc);
1120#else
1121 acc = _mm256_add_ps(_mm256_mul_ps(d, q), acc);
1122#endif
1123 }
1124
1125 *s = hsum_float_8_avx(acc);
1126}
1127#endif
1128
1129#if defined(__SSSE3__)
1130/**
1131 * @brief Spread 32 bits to 32 bytes { 0x00, 0xFF }
1132 * Adapted from llama.cpp bytes_from_bits_32 (AVX path)
1133 *
1134 * Uses shuffle to replicate each byte, then OR with bit_mask and compare.
1135 * Result: 0xFF where bit was set, 0x00 where bit was not set.
1136 */
1137static inline void bytes_from_bits_32_sse(__m128i *out_lo, __m128i *out_hi, const uint8_t *qh)
1138{
1139 uint32_t x32;
1140 memcpy(&x32, qh, sizeof(uint32_t));
1141
1142 /* Shuffle masks: replicate byte j/8 of x32 to each position */
1143 const __m128i shuf_maskl = _mm_set_epi64x(0x0101010101010101LL, 0x0000000000000000LL);
1144 const __m128i shuf_maskh = _mm_set_epi64x(0x0303030303030303LL, 0x0202020202020202LL);
1145
1146 __m128i bytes_lo = _mm_shuffle_epi8(_mm_set1_epi32(x32), shuf_maskl);
1147 __m128i bytes_hi = _mm_shuffle_epi8(_mm_set1_epi32(x32), shuf_maskh);
1148
1149 /* Bit mask: pattern tests each bit position 0-7 within each byte.
1150 * 0x7fbfdfeff7fbfdfe in binary has bits 1,2,3,4,5,6,7,0 cleared per 8-byte cycle.
1151 * After OR, byte will be 0xFF if the corresponding bit was set. */
1152 const __m128i bit_mask = _mm_set1_epi64x(0x7fbfdfeff7fbfdfeLL);
1153
1154 bytes_lo = _mm_or_si128(bytes_lo, bit_mask);
1155 bytes_hi = _mm_or_si128(bytes_hi, bit_mask);
1156
1157 /* Compare with all 1s: 0xFF if bit was set, 0x00 if not */
1158 *out_lo = _mm_cmpeq_epi8(bytes_lo, _mm_set1_epi64x(-1LL));
1159 *out_hi = _mm_cmpeq_epi8(bytes_hi, _mm_set1_epi64x(-1LL));
1160}
1161
1162/**
1163 * @brief Multiply signed int8 vectors using sign trick
1164 * Adapted from llama.cpp mul_sum_i8_pairs
1165 *
1166 * Uses: abs(x) * sign(y,x) = x * y for signed multiplication with maddubs
1167 */
1168static inline __m128i mul_sum_i8_pairs_sse(const __m128i x, const __m128i y)
1169{
1170 const __m128i ax = _mm_sign_epi8(x, x); /* abs(x) */
1171 const __m128i sy = _mm_sign_epi8(y, x); /* y * sign(x) */
1172 const __m128i dot = _mm_maddubs_epi16(ax, sy); /* unsigned*signed pairs -> int16 */
1173 return _mm_madd_epi16(dot, _mm_set1_epi16(1)); /* sum pairs -> int32 */
1174}
1175
1176/**
1177 * @brief Vectorized dot product Q5_0 x Q8_0 using SSSE3
1178 *
1179 * Based on llama.cpp ggml_vec_dot_q5_0_q8_0 AVX implementation.
1180 * Key insight: use shuffle-based bit spreading and sign trick.
1181 *
1182 * Q5_0 encoding: nibble | (high_bit ? 0 : 0xF0)
1183 * - When high bit SET: value = nibble (0-15, positive as signed)
1184 * - When high bit NOT SET: value = nibble | 0xF0 (negative as signed, -16 to -1)
1185 *
1186 * Sign trick handles signed*signed multiplication with unsigned*signed maddubs.
1187 */
1188void vec_dot_q5_0_q8_0_sse(int n, float *s, const void *vx, const void *vy)
1189{
1190 const int qk = QK5_0; /* 32 */
1191 const int nb = n / qk;
1192
1193 const block_q5_0 *x = (const block_q5_0 *)vx;
1194 const block_q8_0 *y = (const block_q8_0 *)vy;
1195
1196 float sumf = 0.0f;
1197
1198 const __m128i mask_0f = _mm_set1_epi8(0x0F);
1199 const __m128i mask_f0 = _mm_set1_epi8((char)0xF0);
1200
1201 for (int ib = 0; ib < nb; ib++) {
1202 const float d = CK_FP16_TO_FP32(x[ib].d) * CK_FP16_TO_FP32(y[ib].d);
1203
1204 /* Load 16 bytes of packed nibbles */
1205 __m128i qs = _mm_loadu_si128((const __m128i *)x[ib].qs);
1206
1207 /* Extract nibbles: lo for indices 0-15, hi for indices 16-31 */
1208 __m128i bx_lo = _mm_and_si128(qs, mask_0f);
1209 __m128i bx_hi = _mm_and_si128(_mm_srli_epi16(qs, 4), mask_0f);
1210
1211 /* Spread 32 high bits to 32 bytes (0xFF=set, 0x00=not set) */
1212 __m128i bxhi_lo, bxhi_hi;
1213 bytes_from_bits_32_sse(&bxhi_lo, &bxhi_hi, x[ib].qh);
1214
1215 /* Apply encoding: (~bxhi) & 0xF0
1216 * When bit SET: bxhi=0xFF, result = 0x00 (value is positive 0-15)
1217 * When bit NOT SET: bxhi=0x00, result = 0xF0 (value is negative) */
1218 bxhi_lo = _mm_andnot_si128(bxhi_lo, mask_f0);
1219 bxhi_hi = _mm_andnot_si128(bxhi_hi, mask_f0);
1220
1221 /* Combine: nibble | high_bit_contribution */
1222 bx_lo = _mm_or_si128(bx_lo, bxhi_lo);
1223 bx_hi = _mm_or_si128(bx_hi, bxhi_hi);
1224
1225 /* Load Q8_0 values (32 signed int8) */
1226 __m128i by_lo = _mm_loadu_si128((const __m128i *)y[ib].qs);
1227 __m128i by_hi = _mm_loadu_si128((const __m128i *)(y[ib].qs + 16));
1228
1229 /* Multiply using sign trick and sum to int32 */
1230 __m128i p_lo = mul_sum_i8_pairs_sse(bx_lo, by_lo);
1231 __m128i p_hi = mul_sum_i8_pairs_sse(bx_hi, by_hi);
1232
1233 /* Sum the two halves */
1234 __m128i sum = _mm_add_epi32(p_lo, p_hi);
1235
1236 /* Horizontal sum of 4 int32 values (avoiding hadd for better perf) */
1237 __m128i hi64 = _mm_unpackhi_epi64(sum, sum);
1238 __m128i sum64 = _mm_add_epi32(hi64, sum);
1239 __m128i hi32 = _mm_shuffle_epi32(sum64, _MM_SHUFFLE(2, 3, 0, 1));
1240 int32_t sumi = _mm_cvtsi128_si32(_mm_add_epi32(sum64, hi32));
1241
1242 /* Scale and accumulate */
1243 sumf += d * (float)sumi;
1244 }
1245
1246 *s = sumf;
1247}
1248#endif
1249
1250#if defined(__AVX__)
1251
1252/* Combine two __m128i into __m256i (AVX without AVX2) */
1253#define MM256_SET_M128I(hi, lo) _mm256_insertf128_si256(_mm256_castsi128_si256(lo), (hi), 1)
1254
1255/**
1256 * @brief Spread 32 bits to 32 bytes using AVX
1257 * Returns __m256i with 0xFF where bit was set, 0x00 where not
1258 */
1259static inline __m256i bytes_from_bits_32_avx(const uint8_t *qh)
1260{
1261 uint32_t x32;
1262 memcpy(&x32, qh, sizeof(uint32_t));
1263
1264 const __m128i shuf_maskl = _mm_set_epi64x(0x0101010101010101LL, 0x0000000000000000LL);
1265 const __m128i shuf_maskh = _mm_set_epi64x(0x0303030303030303LL, 0x0202020202020202LL);
1266
1267 __m128i bytesl = _mm_shuffle_epi8(_mm_set1_epi32(x32), shuf_maskl);
1268 __m128i bytesh = _mm_shuffle_epi8(_mm_set1_epi32(x32), shuf_maskh);
1269
1270 const __m128i bit_mask = _mm_set1_epi64x(0x7fbfdfeff7fbfdfeLL);
1271
1272 bytesl = _mm_or_si128(bytesl, bit_mask);
1273 bytesh = _mm_or_si128(bytesh, bit_mask);
1274
1275 bytesl = _mm_cmpeq_epi8(bytesl, _mm_set1_epi64x(-1LL));
1276 bytesh = _mm_cmpeq_epi8(bytesh, _mm_set1_epi64x(-1LL));
1277
1278 return MM256_SET_M128I(bytesh, bytesl);
1279}
1280
1281/**
1282 * @brief Unpack 32 4-bit nibbles to 32 bytes using AVX
1283 */
1284static inline __m256i bytes_from_nibbles_32_avx(const uint8_t *qs)
1285{
1286 __m128i tmpl = _mm_loadu_si128((const __m128i *)qs);
1287 __m128i tmph = _mm_srli_epi16(tmpl, 4);
1288 const __m128i lowMask = _mm_set1_epi8(0x0F);
1289 tmpl = _mm_and_si128(lowMask, tmpl);
1290 tmph = _mm_and_si128(lowMask, tmph);
1291 return MM256_SET_M128I(tmph, tmpl);
1292}
1293
1294/**
1295 * @brief Multiply signed int8 pairs and return as float vector (AVX)
1296 * Uses 128-bit ops internally but returns 256-bit float result
1297 */
1298static inline __m256 mul_sum_i8_pairs_float_avx(const __m256i x, const __m256i y)
1299{
1300 const __m128i xl = _mm256_castsi256_si128(x);
1301 const __m128i xh = _mm256_extractf128_si256(x, 1);
1302 const __m128i yl = _mm256_castsi256_si128(y);
1303 const __m128i yh = _mm256_extractf128_si256(y, 1);
1304
1305 /* Get absolute values of x vectors */
1306 const __m128i axl = _mm_sign_epi8(xl, xl);
1307 const __m128i axh = _mm_sign_epi8(xh, xh);
1308 /* Sign the values of the y vectors */
1309 const __m128i syl = _mm_sign_epi8(yl, xl);
1310 const __m128i syh = _mm_sign_epi8(yh, xh);
1311
1312 /* Perform multiplication and create 16-bit values */
1313 const __m128i dotl = _mm_maddubs_epi16(axl, syl);
1314 const __m128i doth = _mm_maddubs_epi16(axh, syh);
1315
1316 /* Sum pairs to int32 */
1317 const __m128i ones = _mm_set1_epi16(1);
1318 const __m128i summed_pairsl = _mm_madd_epi16(ones, dotl);
1319 const __m128i summed_pairsh = _mm_madd_epi16(ones, doth);
1320
1321 /* Convert to float */
1322 const __m256i summed_pairs = MM256_SET_M128I(summed_pairsh, summed_pairsl);
1323 return _mm256_cvtepi32_ps(summed_pairs);
1324}
1325
1326/**
1327 * @brief Horizontal sum of 8 floats (AVX)
1328 */
1329static inline float hsum_float_8_avx(const __m256 x)
1330{
1331 __m128 res = _mm256_extractf128_ps(x, 1);
1332 res = _mm_add_ps(res, _mm256_castps256_ps128(x));
1333 res = _mm_add_ps(res, _mm_movehl_ps(res, res));
1334 res = _mm_add_ss(res, _mm_movehdup_ps(res));
1335 return _mm_cvtss_f32(res);
1336}
1337
1338/**
1339 * @brief Quantized dot product Q5_0 x Q8_0 (AVX) - Optimized with 2x unroll
1340 *
1341 * Based on llama.cpp ggml_vec_dot_q5_0_q8_0 AVX implementation.
1342 * Uses 256-bit accumulation and processes 32 values per block.
1343 *
1344 * Optimizations:
1345 * - 2x loop unrolling to reduce loop overhead
1346 * - Prefetching next blocks to hide memory latency
1347 * - Interleaved operations for better instruction-level parallelism
1348 */
1349void vec_dot_q5_0_q8_0_avx(int n, float *s, const void *vx, const void *vy)
1350{
1351 const int qk = QK5_0; /* 32 */
1352 const int nb = n / qk;
1353
1354 const block_q5_0 *x = (const block_q5_0 *)vx;
1355 const block_q8_0 *y = (const block_q8_0 *)vy;
1356
1357 __m256 acc0 = _mm256_setzero_ps();
1358 __m256 acc1 = _mm256_setzero_ps();
1359 const __m128i mask = _mm_set1_epi8((char)0xF0);
1360
1361 /* Process 2 blocks per iteration */
1362 int ib = 0;
1363 for (; ib + 1 < nb; ib += 2) {
1364 /* Prefetch next blocks (2 cache lines ahead) */
1365 _mm_prefetch((const char *)&x[ib + 4], _MM_HINT_T0);
1366 _mm_prefetch((const char *)&y[ib + 4], _MM_HINT_T0);
1367
1368 /* === Block 0 === */
1369 const __m256 d0 = _mm256_set1_ps(CK_FP16_TO_FP32(x[ib].d) * CK_FP16_TO_FP32(y[ib].d));
1370
1371 /* Unpack nibbles to 32 bytes */
1372 __m256i bx0 = bytes_from_nibbles_32_avx(x[ib].qs);
1373
1374 /* Spread high bits */
1375 const __m256i bxhi0 = bytes_from_bits_32_avx(x[ib].qh);
1376 __m128i bxhil0 = _mm256_castsi256_si128(bxhi0);
1377 __m128i bxhih0 = _mm256_extractf128_si256(bxhi0, 1);
1378
1379 /* === Block 1 === (start while block 0 is in flight) */
1380 const __m256 d1 = _mm256_set1_ps(CK_FP16_TO_FP32(x[ib+1].d) * CK_FP16_TO_FP32(y[ib+1].d));
1381
1382 __m256i bx1 = bytes_from_nibbles_32_avx(x[ib+1].qs);
1383 const __m256i bxhi1 = bytes_from_bits_32_avx(x[ib+1].qh);
1384 __m128i bxhil1 = _mm256_castsi256_si128(bxhi1);
1385 __m128i bxhih1 = _mm256_extractf128_si256(bxhi1, 1);
1386
1387 /* === Finish Block 0 === */
1388 bxhil0 = _mm_andnot_si128(bxhil0, mask);
1389 bxhih0 = _mm_andnot_si128(bxhih0, mask);
1390
1391 __m128i bxl0 = _mm256_castsi256_si128(bx0);
1392 __m128i bxh0 = _mm256_extractf128_si256(bx0, 1);
1393 bxl0 = _mm_or_si128(bxl0, bxhil0);
1394 bxh0 = _mm_or_si128(bxh0, bxhih0);
1395 bx0 = MM256_SET_M128I(bxh0, bxl0);
1396
1397 const __m256i by0 = _mm256_loadu_si256((const __m256i *)y[ib].qs);
1398 const __m256 q0 = mul_sum_i8_pairs_float_avx(bx0, by0);
1399 acc0 = _mm256_add_ps(_mm256_mul_ps(d0, q0), acc0);
1400
1401 /* === Finish Block 1 === */
1402 bxhil1 = _mm_andnot_si128(bxhil1, mask);
1403 bxhih1 = _mm_andnot_si128(bxhih1, mask);
1404
1405 __m128i bxl1 = _mm256_castsi256_si128(bx1);
1406 __m128i bxh1 = _mm256_extractf128_si256(bx1, 1);
1407 bxl1 = _mm_or_si128(bxl1, bxhil1);
1408 bxh1 = _mm_or_si128(bxh1, bxhih1);
1409 bx1 = MM256_SET_M128I(bxh1, bxl1);
1410
1411 const __m256i by1 = _mm256_loadu_si256((const __m256i *)y[ib+1].qs);
1412 const __m256 q1 = mul_sum_i8_pairs_float_avx(bx1, by1);
1413 acc1 = _mm256_add_ps(_mm256_mul_ps(d1, q1), acc1);
1414 }
1415
1416 /* Handle remaining block if nb is odd */
1417 for (; ib < nb; ib++) {
1418 const __m256 d = _mm256_set1_ps(CK_FP16_TO_FP32(x[ib].d) * CK_FP16_TO_FP32(y[ib].d));
1419
1420 __m256i bx_0 = bytes_from_nibbles_32_avx(x[ib].qs);
1421 const __m256i bxhi = bytes_from_bits_32_avx(x[ib].qh);
1422 __m128i bxhil = _mm256_castsi256_si128(bxhi);
1423 __m128i bxhih = _mm256_extractf128_si256(bxhi, 1);
1424
1425 bxhil = _mm_andnot_si128(bxhil, mask);
1426 bxhih = _mm_andnot_si128(bxhih, mask);
1427
1428 __m128i bxl = _mm256_castsi256_si128(bx_0);
1429 __m128i bxh = _mm256_extractf128_si256(bx_0, 1);
1430 bxl = _mm_or_si128(bxl, bxhil);
1431 bxh = _mm_or_si128(bxh, bxhih);
1432 bx_0 = MM256_SET_M128I(bxh, bxl);
1433
1434 const __m256i by_0 = _mm256_loadu_si256((const __m256i *)y[ib].qs);
1435 const __m256 q = mul_sum_i8_pairs_float_avx(bx_0, by_0);
1436 acc0 = _mm256_add_ps(_mm256_mul_ps(d, q), acc0);
1437 }
1438
1439 /* Combine accumulators and sum */
1440 acc0 = _mm256_add_ps(acc0, acc1);
1441 *s = hsum_float_8_avx(acc0);
1442}
1443
1444/* ============================================================================
1445 * Unrolled GEMM: C[M,N] = A_q8[M,K] @ B_q5[N,K]^T + bias
1446 *
1447 * Key optimizations over the naive vec_dot-per-element approach:
1448 * 1. N-dimension 2x unroll: loads activation block ONCE, reuses for 2 weight rows
1449 * 2. Inlined Q5_0 unpack: eliminates function call overhead per block
1450 * 3. Block-level prefetching: hides DRAM latency for weight streaming
1451 * 4. F16C scale conversion: uses hardware vcvtsh2ss via CK_FP16_TO_FP32 macro
1452 * 5. Persistent __m256 accumulators: no per-block horizontal sum
1453 *
1454 * Expected speedup: 1.5-2x over naive dispatch loop on AVX (Ivy Bridge+)
1455 * ============================================================================ */
1456
1457/**
1458 * @brief Batch GEMM with Q5_0 weights x Q8_0 activations — AVX unrolled
1459 *
1460 * @param A_q8 Input activations in Q8_0 format [M rows of K/32 blocks each]
1461 * @param B_q5 Weights in Q5_0 format [N rows of K/32 blocks each]
1462 * @param bias Optional bias vector [N], NULL if not used
1463 * @param C Output matrix [M x N], row-major FP32
1464 * @param M Batch size (number of tokens)
1465 * @param N Output dimension (number of output features)
1466 * @param K Input dimension (must be multiple of 32)
1467 */
1468void gemm_nt_q5_0_q8_0_unroll_avx(
1469 const void *A_q8,
1470 const void *B_q5,
1471 const float *bias,
1472 float *C,
1473 int M, int N, int K)
1474{
1475 const int nb = K / QK5_0;
1476 const block_q8_0 *a_blocks = (const block_q8_0 *)A_q8;
1477 const block_q5_0 *b_blocks = (const block_q5_0 *)B_q5;
1478 const __m128i mask = _mm_set1_epi8((char)0xF0);
1479
1480 for (int m = 0; m < M; m++) {
1481 const block_q8_0 *a_row = a_blocks + (size_t)m * nb;
1482
1483 /* ---- 2x N-unrolled path: process 2 weight rows per iteration ---- */
1484 int n = 0;
1485 for (; n + 1 < N; n += 2) {
1486 const block_q5_0 *w0 = b_blocks + (size_t)(n + 0) * nb;
1487 const block_q5_0 *w1 = b_blocks + (size_t)(n + 1) * nb;
1488
1489 __m256 acc_n0 = _mm256_setzero_ps();
1490 __m256 acc_n1 = _mm256_setzero_ps();
1491
1492 for (int ib = 0; ib < nb; ib++) {
1493 /* Prefetch next weight blocks (1 cache line = ~2.9 Q5_0 blocks) */
1494 if (ib + 2 < nb) {
1495 _mm_prefetch((const char *)&w0[ib + 2], _MM_HINT_T0);
1496 _mm_prefetch((const char *)&w1[ib + 2], _MM_HINT_T0);
1497 }
1498
1499 /* Load activation Q8_0 block ONCE — reused for both weight rows */
1500 const __m256i by = _mm256_loadu_si256((const __m256i *)a_row[ib].qs);
1501 const float da = CK_FP16_TO_FP32(a_row[ib].d);
1502
1503 /* === Weight row 0: inline Q5_0 unpack + dot product === */
1504 {
1505 const float d = da * CK_FP16_TO_FP32(w0[ib].d);
1506
1507 /* Unpack low 4-bit nibbles → 32 bytes */
1508 __m256i bx = bytes_from_nibbles_32_avx(w0[ib].qs);
1509
1510 /* Spread 32 high bits → 32 byte mask (0xFF or 0x00) */
1511 const __m256i bxhi = bytes_from_bits_32_avx(w0[ib].qh);
1512 __m128i bxhil = _mm256_castsi256_si128(bxhi);
1513 __m128i bxhih = _mm256_extractf128_si256(bxhi, 1);
1514
1515 /* Combine: nibble |= (highbit ? 0x10 : 0x00) */
1516 bxhil = _mm_andnot_si128(bxhil, mask);
1517 bxhih = _mm_andnot_si128(bxhih, mask);
1518 __m128i bxl = _mm256_castsi256_si128(bx);
1519 __m128i bxh = _mm256_extractf128_si256(bx, 1);
1520 bxl = _mm_or_si128(bxl, bxhil);
1521 bxh = _mm_or_si128(bxh, bxhih);
1522 bx = MM256_SET_M128I(bxh, bxl);
1523
1524 /* Signed int8 dot product → 8 float partial sums */
1525 const __m256 q = mul_sum_i8_pairs_float_avx(bx, by);
1526 acc_n0 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(d), q), acc_n0);
1527 }
1528
1529 /* === Weight row 1: same activation data, different weights === */
1530 {
1531 const float d = da * CK_FP16_TO_FP32(w1[ib].d);
1532
1533 __m256i bx = bytes_from_nibbles_32_avx(w1[ib].qs);
1534 const __m256i bxhi = bytes_from_bits_32_avx(w1[ib].qh);
1535 __m128i bxhil = _mm256_castsi256_si128(bxhi);
1536 __m128i bxhih = _mm256_extractf128_si256(bxhi, 1);
1537
1538 bxhil = _mm_andnot_si128(bxhil, mask);
1539 bxhih = _mm_andnot_si128(bxhih, mask);
1540 __m128i bxl = _mm256_castsi256_si128(bx);
1541 __m128i bxh = _mm256_extractf128_si256(bx, 1);
1542 bxl = _mm_or_si128(bxl, bxhil);
1543 bxh = _mm_or_si128(bxh, bxhih);
1544 bx = MM256_SET_M128I(bxh, bxl);
1545
1546 const __m256 q = mul_sum_i8_pairs_float_avx(bx, by);
1547 acc_n1 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(d), q), acc_n1);
1548 }
1549 }
1550
1551 /* Horizontal reduce and store with optional bias */
1552 float s0 = hsum_float_8_avx(acc_n0);
1553 float s1 = hsum_float_8_avx(acc_n1);
1554 if (bias) { s0 += bias[n]; s1 += bias[n + 1]; }
1555 C[(size_t)m * N + n] = s0;
1556 C[(size_t)m * N + n + 1] = s1;
1557 }
1558
1559 /* ---- Cleanup: remaining odd column ---- */
1560 for (; n < N; n++) {
1561 const block_q5_0 *w = b_blocks + (size_t)n * nb;
1562 __m256 acc = _mm256_setzero_ps();
1563
1564 for (int ib = 0; ib < nb; ib++) {
1565 const __m256i by = _mm256_loadu_si256((const __m256i *)a_row[ib].qs);
1566 const float d = CK_FP16_TO_FP32(a_row[ib].d) * CK_FP16_TO_FP32(w[ib].d);
1567
1568 __m256i bx = bytes_from_nibbles_32_avx(w[ib].qs);
1569 const __m256i bxhi = bytes_from_bits_32_avx(w[ib].qh);
1570 __m128i bxhil = _mm256_castsi256_si128(bxhi);
1571 __m128i bxhih = _mm256_extractf128_si256(bxhi, 1);
1572 bxhil = _mm_andnot_si128(bxhil, mask);
1573 bxhih = _mm_andnot_si128(bxhih, mask);
1574 __m128i bxl = _mm256_castsi256_si128(bx);
1575 __m128i bxh = _mm256_extractf128_si256(bx, 1);
1576 bxl = _mm_or_si128(bxl, bxhil);
1577 bxh = _mm_or_si128(bxh, bxhih);
1578 bx = MM256_SET_M128I(bxh, bxl);
1579
1580 const __m256 q = mul_sum_i8_pairs_float_avx(bx, by);
1581 acc = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(d), q), acc);
1582 }
1583
1584 float s = hsum_float_8_avx(acc);
1585 if (bias) s += bias[n];
1586 C[(size_t)m * N + n] = s;
1587 }
1588 }
1589}
1590
1591#endif
1592
1593/**
1594 * @brief Auto-dispatch quantized dot product Q5_0 x Q8_0
1595 *
1596 * Dispatch priority:
1597 * 1. AVX512 (best performance on modern Intel/AMD)
1598 * 2. AVX (256-bit float ops, works on Sandy/Ivy Bridge and newer)
1599 * 3. SSSE3 (128-bit fallback)
1600 * 4. Reference scalar (last resort)
1601 */
1602void vec_dot_q5_0_q8_0(int n, float *s, const void *vx, const void *vy)
1603{
1604#if defined(__AVX2__)
1605 /* llama.cpp uses the packed AVX2 dot on AVX-512 hosts as well. It keeps
1606 * Q5/Q8 data in byte lanes and avoids the per-block 32-bit lane expansion
1607 * overhead of the baseline AVX-512 path. */
1608 vec_dot_q5_0_q8_0_avx2(n, s, vx, vy);
1609#elif defined(__AVX512F__)
1610 vec_dot_q5_0_q8_0_avx512(n, s, vx, vy);
1611#elif defined(__ARM_NEON) || defined(__aarch64__)
1612 vec_dot_q5_0_q8_0_neon(n, s, vx, vy);
1613#elif defined(__AVX__)
1614 /* AVX for 256-bit float ops (works on Ivy Bridge and newer) */
1615 vec_dot_q5_0_q8_0_avx(n, s, vx, vy);
1616#elif defined(__SSSE3__)
1617 /* SSSE3 - most efficient on older CPUs */
1618 vec_dot_q5_0_q8_0_sse(n, s, vx, vy);
1619#else
1620 vec_dot_q5_0_q8_0_ref(n, s, vx, vy);
1621#endif
1622}
1623
1624/* ============================================================================
1625 * Quantized GEMV: y = W @ x where W is Q5_0 and x is Q8_0
1626 *
1627 * This is the quantized equivalent of gemv_q5_0, but takes pre-quantized
1628 * input in Q8_0 format. Used for parity testing with llama.cpp.
1629 * ============================================================================ */
1630
1631/**
1632 * @brief Matrix-vector multiply with Q5_0 weights and Q8_0 input
1633 *
1634 * @param y Output vector [M]
1635 * @param W Weight matrix in Q5_0 format [M x K]
1636 * @param x_q8 Input vector in Q8_0 format [K]
1637 * @param M Number of output rows
1638 * @param K Number of columns (must be multiple of 32)
1639 */
1640void gemv_q5_0_q8_0(float *y,
1641 const void *W,
1642 const void *x_q8,
1643 int M, int K)
1644{
1645 const block_q5_0 *w_blocks = (const block_q5_0 *)W;
1646 const block_q8_0 *x_blocks = (const block_q8_0 *)x_q8;
1647 const int blocks_per_row = K / QK5_0;
1648
1649 for (int row = 0; row < M; row++) {
1650 vec_dot_q5_0_q8_0(K, &y[row],
1651 &w_blocks[row * blocks_per_row],
1652 x_blocks);
1653 }
1654}
1655
1656/**
1657 * @brief Parallel SIMD GEMV for Q5_0 x Q8_0 with prefetching
1658 *
1659 * Each thread processes rows [r0, r1) where r0 = ith * ceil(M/nth).
1660 * Uses vec_dot_q5_0_q8_0 dispatch (auto-selects AVX512/AVX/SSE/scalar).
1661 */
1662void gemv_q5_0_q8_0_parallel_simd(float *y,
1663 const void *W,
1664 const void *x_q8,
1665 int M, int K,
1666 int ith, int nth)
1667{
1668 if (!y || !W || !x_q8 || M <= 0 || K <= 0) return;
1669 if (ith < 0 || nth <= 0 || ith >= nth) return;
1670
1671 const int dr = (M + nth - 1) / nth;
1672 const int r0 = dr * ith;
1673 const int r1 = (r0 + dr < M) ? (r0 + dr) : M;
1674
1675 if (r0 >= M) return;
1676
1677 const block_q5_0 *w_blocks = (const block_q5_0 *)W;
1678 const block_q8_0 *x_blocks = (const block_q8_0 *)x_q8;
1679 const int blocks_per_row = K / QK5_0;
1680
1681#if defined(__AVX__) || defined(__SSE4_1__)
1682 const int PREFETCH_ROWS = 4;
1683 for (int p = 0; p < PREFETCH_ROWS && r0 + p < r1; ++p) {
1684 const char *row_ptr = (const char *)(w_blocks + (r0 + p) * blocks_per_row);
1685 _mm_prefetch(row_ptr, _MM_HINT_T0);
1686 _mm_prefetch(row_ptr + 64, _MM_HINT_T0);
1687 }
1688
1689 for (int row = r0; row < r1; ++row) {
1690 if (row + PREFETCH_ROWS < r1) {
1691 const char *pf = (const char *)(w_blocks + (row + PREFETCH_ROWS) * blocks_per_row);
1692 _mm_prefetch(pf, _MM_HINT_T0);
1693 _mm_prefetch(pf + 64, _MM_HINT_T0);
1694 }
1695
1696 vec_dot_q5_0_q8_0(K, &y[row],
1697 &w_blocks[row * blocks_per_row],
1698 x_blocks);
1699 }
1700#else
1701 for (int row = r0; row < r1; row++) {
1702 vec_dot_q5_0_q8_0(K, &y[row],
1703 &w_blocks[row * blocks_per_row],
1704 x_blocks);
1705 }
1706#endif
1707}
1708
1709/**
1710 * @brief Batch GEMM with Q5_0 weights and Q8_0 activations for prefill
1711 *
1712 * Computes C = A @ B^T + bias where:
1713 * A: [M x K] Q8_0 quantized activations (M tokens, K features)
1714 * B: [N x K] Q5_0 quantized weights (N outputs, K features)
1715 * C: [M x N] FP32 output
1716 *
1717 * This is the INT8 batch kernel for prefill, using pre-quantized activations
1718 * to avoid FP32->Q8_0 conversion overhead per operation.
1719 *
1720 * @param A_q8 Input activations in Q8_0 format [M rows of K/32 blocks each]
1721 * @param B_q5 Weights in Q5_0 format [N rows of K/32 blocks each]
1722 * @param bias Optional bias vector [N], NULL if not used
1723 * @param C Output matrix [M x N], row-major FP32
1724 * @param M Batch size (number of tokens)
1725 * @param N Output dimension (number of output features)
1726 * @param K Input dimension (must be multiple of 32)
1727 */
1728void gemm_nt_q5_0_q8_0(
1729 const void *A_q8,
1730 const void *B_q5,
1731 const float *bias,
1732 float *C,
1733 int M,
1734 int N,
1735 int K)
1736{
1737 const block_q5_0 *weights = (const block_q5_0 *)B_q5;
1738 const block_q8_0 *inputs = (const block_q8_0 *)A_q8;
1739 const int blocks_per_row = K / QK5_0;
1740
1741 for (int m = 0; m < M; m++) {
1742 const block_q8_0 *input_row = &inputs[m * blocks_per_row];
1743
1744 for (int n = 0; n < N; n++) {
1745 const block_q5_0 *weight_row = &weights[n * blocks_per_row];
1746 float *out = &C[m * N + n];
1747
1748 /* Dispatches to vec_dot_q5_0_q8_0_avx (2x block unrolled) on AVX */
1749 vec_dot_q5_0_q8_0(K, out, weight_row, input_row);
1750
1751 if (bias) {
1752 *out += bias[n];
1753 }
1754 }
1755 }
1756}
1757
1758/*
1759 * Two-token by four-output Q5_0 provider.
1760 *
1761 * Q5 unpacking is independent of the activation row. Keep the certified
1762 * eight-lane accumulator and horizontal reduction for every output while
1763 * applying each unpacked weight block to two token rows. This removes one
1764 * complete Q5 reconstruction per token pair without coupling reductions.
1765 */
1766void gemm_nt_q5_0_q8_0_m2n4_tile(
1767 const void *A_q8,
1768 const void *B_q5,
1769 const float *bias,
1770 float *C,
1771 int M,
1772 int N,
1773 int K,
1774 int ldc)
1775{
1776#if defined(__AVX2__)
1777 if (!A_q8 || !B_q5 || !C || M <= 0 || N <= 0 || K <= 0 ||
1778 (K % QK5_0) != 0 || ldc < N) {
1779 return;
1780 }
1781
1782 const block_q8_0 *a = (const block_q8_0 *)A_q8;
1783 const block_q5_0 *w = (const block_q5_0 *)B_q5;
1784 const int nb = K / QK5_0;
1785 int m = 0;
1786
1787 for (; m + 1 < M; m += 2) {
1788 const block_q8_0 *a0 = a + (size_t)(m + 0) * (size_t)nb;
1789 const block_q8_0 *a1 = a + (size_t)(m + 1) * (size_t)nb;
1790 int n = 0;
1791
1792 for (; n + 3 < N; n += 4) {
1793 const block_q5_0 *w0 = w + (size_t)(n + 0) * (size_t)nb;
1794 const block_q5_0 *w1 = w + (size_t)(n + 1) * (size_t)nb;
1795 const block_q5_0 *w2 = w + (size_t)(n + 2) * (size_t)nb;
1796 const block_q5_0 *w3 = w + (size_t)(n + 3) * (size_t)nb;
1797 __m256 acc00 = _mm256_setzero_ps();
1798 __m256 acc01 = _mm256_setzero_ps();
1799 __m256 acc02 = _mm256_setzero_ps();
1800 __m256 acc03 = _mm256_setzero_ps();
1801 __m256 acc10 = _mm256_setzero_ps();
1802 __m256 acc11 = _mm256_setzero_ps();
1803 __m256 acc12 = _mm256_setzero_ps();
1804 __m256 acc13 = _mm256_setzero_ps();
1805
1806 for (int ib = 0; ib < nb; ++ib) {
1807 __m256i qw0 = bytes_from_nibbles_32_avx(w0[ib].qs);
1808 __m256i qw1 = bytes_from_nibbles_32_avx(w1[ib].qs);
1809 __m256i qw2 = bytes_from_nibbles_32_avx(w2[ib].qs);
1810 __m256i qw3 = bytes_from_nibbles_32_avx(w3[ib].qs);
1811 const __m256i sign = _mm256_set1_epi8((char)0xF0);
1812 qw0 = _mm256_or_si256(qw0, _mm256_andnot_si256(bytes_from_bits_32_avx(w0[ib].qh), sign));
1813 qw1 = _mm256_or_si256(qw1, _mm256_andnot_si256(bytes_from_bits_32_avx(w1[ib].qh), sign));
1814 qw2 = _mm256_or_si256(qw2, _mm256_andnot_si256(bytes_from_bits_32_avx(w2[ib].qh), sign));
1815 qw3 = _mm256_or_si256(qw3, _mm256_andnot_si256(bytes_from_bits_32_avx(w3[ib].qh), sign));
1816
1817 const __m256i qa0 = _mm256_loadu_si256((const __m256i *)a0[ib].qs);
1818 const __m256i qa1 = _mm256_loadu_si256((const __m256i *)a1[ib].qs);
1819 const float da0 = CK_FP16_TO_FP32(a0[ib].d);
1820 const float da1 = CK_FP16_TO_FP32(a1[ib].d);
1821 const float dw0 = CK_FP16_TO_FP32(w0[ib].d);
1822 const float dw1 = CK_FP16_TO_FP32(w1[ib].d);
1823 const float dw2 = CK_FP16_TO_FP32(w2[ib].d);
1824 const float dw3 = CK_FP16_TO_FP32(w3[ib].d);
1825 const __m256 p00 = mul_sum_i8_pairs_float_avx(qw0, qa0);
1826 const __m256 p01 = mul_sum_i8_pairs_float_avx(qw1, qa0);
1827 const __m256 p02 = mul_sum_i8_pairs_float_avx(qw2, qa0);
1828 const __m256 p03 = mul_sum_i8_pairs_float_avx(qw3, qa0);
1829 const __m256 p10 = mul_sum_i8_pairs_float_avx(qw0, qa1);
1830 const __m256 p11 = mul_sum_i8_pairs_float_avx(qw1, qa1);
1831 const __m256 p12 = mul_sum_i8_pairs_float_avx(qw2, qa1);
1832 const __m256 p13 = mul_sum_i8_pairs_float_avx(qw3, qa1);
1833#if defined(__FMA__)
1834 acc00 = _mm256_fmadd_ps(_mm256_set1_ps(dw0 * da0), p00, acc00);
1835 acc01 = _mm256_fmadd_ps(_mm256_set1_ps(dw1 * da0), p01, acc01);
1836 acc02 = _mm256_fmadd_ps(_mm256_set1_ps(dw2 * da0), p02, acc02);
1837 acc03 = _mm256_fmadd_ps(_mm256_set1_ps(dw3 * da0), p03, acc03);
1838 acc10 = _mm256_fmadd_ps(_mm256_set1_ps(dw0 * da1), p10, acc10);
1839 acc11 = _mm256_fmadd_ps(_mm256_set1_ps(dw1 * da1), p11, acc11);
1840 acc12 = _mm256_fmadd_ps(_mm256_set1_ps(dw2 * da1), p12, acc12);
1841 acc13 = _mm256_fmadd_ps(_mm256_set1_ps(dw3 * da1), p13, acc13);
1842#else
1843 acc00 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw0 * da0), p00), acc00);
1844 acc01 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw1 * da0), p01), acc01);
1845 acc02 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw2 * da0), p02), acc02);
1846 acc03 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw3 * da0), p03), acc03);
1847 acc10 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw0 * da1), p10), acc10);
1848 acc11 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw1 * da1), p11), acc11);
1849 acc12 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw2 * da1), p12), acc12);
1850 acc13 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw3 * da1), p13), acc13);
1851#endif
1852 }
1853
1854 C[(size_t)(m + 0) * (size_t)ldc + n + 0] = hsum_float_8_avx(acc00) + (bias ? bias[n + 0] : 0.0f);
1855 C[(size_t)(m + 0) * (size_t)ldc + n + 1] = hsum_float_8_avx(acc01) + (bias ? bias[n + 1] : 0.0f);
1856 C[(size_t)(m + 0) * (size_t)ldc + n + 2] = hsum_float_8_avx(acc02) + (bias ? bias[n + 2] : 0.0f);
1857 C[(size_t)(m + 0) * (size_t)ldc + n + 3] = hsum_float_8_avx(acc03) + (bias ? bias[n + 3] : 0.0f);
1858 C[(size_t)(m + 1) * (size_t)ldc + n + 0] = hsum_float_8_avx(acc10) + (bias ? bias[n + 0] : 0.0f);
1859 C[(size_t)(m + 1) * (size_t)ldc + n + 1] = hsum_float_8_avx(acc11) + (bias ? bias[n + 1] : 0.0f);
1860 C[(size_t)(m + 1) * (size_t)ldc + n + 2] = hsum_float_8_avx(acc12) + (bias ? bias[n + 2] : 0.0f);
1861 C[(size_t)(m + 1) * (size_t)ldc + n + 3] = hsum_float_8_avx(acc13) + (bias ? bias[n + 3] : 0.0f);
1862 }
1863
1864 for (; n < N; ++n) {
1865 const block_q5_0 *wn = w + (size_t)n * (size_t)nb;
1866 vec_dot_q5_0_q8_0(K, &C[(size_t)(m + 0) * (size_t)ldc + n], wn, a0);
1867 vec_dot_q5_0_q8_0(K, &C[(size_t)(m + 1) * (size_t)ldc + n], wn, a1);
1868 if (bias) {
1869 C[(size_t)(m + 0) * (size_t)ldc + n] += bias[n];
1870 C[(size_t)(m + 1) * (size_t)ldc + n] += bias[n];
1871 }
1872 }
1873 }
1874
1875 if (m < M) {
1876 gemm_nt_q5_0_q8_0(
1877 a + (size_t)m * (size_t)nb, w, bias,
1878 C + (size_t)m * (size_t)ldc, 1, N, K);
1879 }
1880#else
1881 if (ldc == N) {
1882 gemm_nt_q5_0_q8_0(A_q8, B_q5, bias, C, M, N, K);
1883 } else {
1884 const int nb = K / QK5_0;
1885 const block_q8_0 *a = (const block_q8_0 *)A_q8;
1886 for (int m = 0; m < M; ++m) {
1887 gemm_nt_q5_0_q8_0(a + (size_t)m * (size_t)nb, B_q5, bias,
1888 C + (size_t)m * (size_t)ldc, 1, N, K);
1889 }
1890 }
1891#endif
1892}
1893
1894void gemm_nt_q5_0_q8_0_m2n4(
1895 const void *A_q8,
1896 const void *B_q5,
1897 const float *bias,
1898 float *C,
1899 int M,
1900 int N,
1901 int K)
1902{
1903 gemm_nt_q5_0_q8_0_m2n4_tile(A_q8, B_q5, bias, C, M, N, K, N);
1904}
1905
1906/* Four-token by two-output companion to m2n4. It keeps the same eight
1907 * independent accumulators but amortizes Q5 reconstruction over four rows. */
1908void gemm_nt_q5_0_q8_0_m4n2_tile(
1909 const void *A_q8,
1910 const void *B_q5,
1911 const float *bias,
1912 float *C,
1913 int M,
1914 int N,
1915 int K,
1916 int ldc)
1917{
1918#if defined(__AVX2__)
1919 if (!A_q8 || !B_q5 || !C || M <= 0 || N <= 0 || K <= 0 ||
1920 (K % QK5_0) != 0 || ldc < N) {
1921 return;
1922 }
1923
1924 const block_q8_0 *a = (const block_q8_0 *)A_q8;
1925 const block_q5_0 *w = (const block_q5_0 *)B_q5;
1926 const int nb = K / QK5_0;
1927 int m = 0;
1928
1929 for (; m + 3 < M; m += 4) {
1930 const block_q8_0 *a0 = a + (size_t)(m + 0) * (size_t)nb;
1931 const block_q8_0 *a1 = a + (size_t)(m + 1) * (size_t)nb;
1932 const block_q8_0 *a2 = a + (size_t)(m + 2) * (size_t)nb;
1933 const block_q8_0 *a3 = a + (size_t)(m + 3) * (size_t)nb;
1934 int n = 0;
1935
1936 for (; n + 1 < N; n += 2) {
1937 const block_q5_0 *w0 = w + (size_t)(n + 0) * (size_t)nb;
1938 const block_q5_0 *w1 = w + (size_t)(n + 1) * (size_t)nb;
1939 __m256 acc00 = _mm256_setzero_ps();
1940 __m256 acc01 = _mm256_setzero_ps();
1941 __m256 acc10 = _mm256_setzero_ps();
1942 __m256 acc11 = _mm256_setzero_ps();
1943 __m256 acc20 = _mm256_setzero_ps();
1944 __m256 acc21 = _mm256_setzero_ps();
1945 __m256 acc30 = _mm256_setzero_ps();
1946 __m256 acc31 = _mm256_setzero_ps();
1947
1948 for (int ib = 0; ib < nb; ++ib) {
1949 __m256i qw0 = bytes_from_nibbles_32_avx(w0[ib].qs);
1950 __m256i qw1 = bytes_from_nibbles_32_avx(w1[ib].qs);
1951 const __m256i sign = _mm256_set1_epi8((char)0xF0);
1952 qw0 = _mm256_or_si256(qw0, _mm256_andnot_si256(bytes_from_bits_32_avx(w0[ib].qh), sign));
1953 qw1 = _mm256_or_si256(qw1, _mm256_andnot_si256(bytes_from_bits_32_avx(w1[ib].qh), sign));
1954
1955 const __m256i qa0 = _mm256_loadu_si256((const __m256i *)a0[ib].qs);
1956 const __m256i qa1 = _mm256_loadu_si256((const __m256i *)a1[ib].qs);
1957 const __m256i qa2 = _mm256_loadu_si256((const __m256i *)a2[ib].qs);
1958 const __m256i qa3 = _mm256_loadu_si256((const __m256i *)a3[ib].qs);
1959 const float dw0 = CK_FP16_TO_FP32(w0[ib].d);
1960 const float dw1 = CK_FP16_TO_FP32(w1[ib].d);
1961 const float da0 = CK_FP16_TO_FP32(a0[ib].d);
1962 const float da1 = CK_FP16_TO_FP32(a1[ib].d);
1963 const float da2 = CK_FP16_TO_FP32(a2[ib].d);
1964 const float da3 = CK_FP16_TO_FP32(a3[ib].d);
1965#if defined(__FMA__)
1966 acc00 = _mm256_fmadd_ps(_mm256_set1_ps(dw0 * da0), mul_sum_i8_pairs_float_avx(qw0, qa0), acc00);
1967 acc01 = _mm256_fmadd_ps(_mm256_set1_ps(dw1 * da0), mul_sum_i8_pairs_float_avx(qw1, qa0), acc01);
1968 acc10 = _mm256_fmadd_ps(_mm256_set1_ps(dw0 * da1), mul_sum_i8_pairs_float_avx(qw0, qa1), acc10);
1969 acc11 = _mm256_fmadd_ps(_mm256_set1_ps(dw1 * da1), mul_sum_i8_pairs_float_avx(qw1, qa1), acc11);
1970 acc20 = _mm256_fmadd_ps(_mm256_set1_ps(dw0 * da2), mul_sum_i8_pairs_float_avx(qw0, qa2), acc20);
1971 acc21 = _mm256_fmadd_ps(_mm256_set1_ps(dw1 * da2), mul_sum_i8_pairs_float_avx(qw1, qa2), acc21);
1972 acc30 = _mm256_fmadd_ps(_mm256_set1_ps(dw0 * da3), mul_sum_i8_pairs_float_avx(qw0, qa3), acc30);
1973 acc31 = _mm256_fmadd_ps(_mm256_set1_ps(dw1 * da3), mul_sum_i8_pairs_float_avx(qw1, qa3), acc31);
1974#else
1975 acc00 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw0 * da0), mul_sum_i8_pairs_float_avx(qw0, qa0)), acc00);
1976 acc01 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw1 * da0), mul_sum_i8_pairs_float_avx(qw1, qa0)), acc01);
1977 acc10 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw0 * da1), mul_sum_i8_pairs_float_avx(qw0, qa1)), acc10);
1978 acc11 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw1 * da1), mul_sum_i8_pairs_float_avx(qw1, qa1)), acc11);
1979 acc20 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw0 * da2), mul_sum_i8_pairs_float_avx(qw0, qa2)), acc20);
1980 acc21 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw1 * da2), mul_sum_i8_pairs_float_avx(qw1, qa2)), acc21);
1981 acc30 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw0 * da3), mul_sum_i8_pairs_float_avx(qw0, qa3)), acc30);
1982 acc31 = _mm256_add_ps(_mm256_mul_ps(_mm256_set1_ps(dw1 * da3), mul_sum_i8_pairs_float_avx(qw1, qa3)), acc31);
1983#endif
1984 }
1985
1986 C[(size_t)(m + 0) * (size_t)ldc + n + 0] = hsum_float_8_avx(acc00) + (bias ? bias[n + 0] : 0.0f);
1987 C[(size_t)(m + 0) * (size_t)ldc + n + 1] = hsum_float_8_avx(acc01) + (bias ? bias[n + 1] : 0.0f);
1988 C[(size_t)(m + 1) * (size_t)ldc + n + 0] = hsum_float_8_avx(acc10) + (bias ? bias[n + 0] : 0.0f);
1989 C[(size_t)(m + 1) * (size_t)ldc + n + 1] = hsum_float_8_avx(acc11) + (bias ? bias[n + 1] : 0.0f);
1990 C[(size_t)(m + 2) * (size_t)ldc + n + 0] = hsum_float_8_avx(acc20) + (bias ? bias[n + 0] : 0.0f);
1991 C[(size_t)(m + 2) * (size_t)ldc + n + 1] = hsum_float_8_avx(acc21) + (bias ? bias[n + 1] : 0.0f);
1992 C[(size_t)(m + 3) * (size_t)ldc + n + 0] = hsum_float_8_avx(acc30) + (bias ? bias[n + 0] : 0.0f);
1993 C[(size_t)(m + 3) * (size_t)ldc + n + 1] = hsum_float_8_avx(acc31) + (bias ? bias[n + 1] : 0.0f);
1994 }
1995
1996 if (n < N) {
1997 const block_q5_0 *wn = w + (size_t)n * (size_t)nb;
1998 const block_q8_0 *rows[4] = {a0, a1, a2, a3};
1999 for (int r = 0; r < 4; ++r) {
2000 float *out = &C[(size_t)(m + r) * (size_t)ldc + n];
2001 vec_dot_q5_0_q8_0(K, out, wn, rows[r]);
2002 if (bias) *out += bias[n];
2003 }
2004 }
2005 }
2006
2007 if (m < M) {
2008 gemm_nt_q5_0_q8_0_m2n4_tile(
2009 a + (size_t)m * (size_t)nb, w, bias,
2010 C + (size_t)m * (size_t)ldc, M - m, N, K, ldc);
2011 }
2012#else
2013 gemm_nt_q5_0_q8_0_m2n4_tile(A_q8, B_q5, bias, C, M, N, K, ldc);
2014#endif
2015}
2016
2017void gemm_nt_q5_0_q8_0_m4n2(
2018 const void *A_q8,
2019 const void *B_q5,
2020 const float *bias,
2021 float *C,
2022 int M,
2023 int N,
2024 int K)
2025{
2026 gemm_nt_q5_0_q8_0_m4n2_tile(A_q8, B_q5, bias, C, M, N, K, N);
2027}
ck_features.h
CPU feature detection and dispatch macros.
ckernel_quant.h
Quantization block structures for weight-only quantization.
QK5_0
#define QK5_0
Definition ckernel_quant.h:68
CK_FP16_TO_FP32
#define CK_FP16_TO_FP32(x)
Definition ckernel_quant.h:440
gemm_nt_q5_0_q8_0_unroll_avx
void gemm_nt_q5_0_q8_0_unroll_avx(const void *A_q8, const void *B_q5, const float *bias, float *C, int M, int N, int K)
gemm_nt_q5_0_sse_v2
void gemm_nt_q5_0_sse_v2(const float *A, const void *B, const float *bias, float *C, int M, int N, int K)
Definition gemm_kernels_q5_0_sse_v2.c:77
gemm_q5_0_backward
void gemm_q5_0_backward(float *dX, const void *W, const float *dY, int M, int N, int K)
Batched backward pass.
Definition gemm_kernels_q5_0.c:764
gemv_q5_0_backward_ref
void gemv_q5_0_backward_ref(float *dX, const void *W, const float *dY, int M, int K)
Backward pass: compute input gradient.
Definition gemm_kernels_q5_0.c:707
dequant_q5_0_block
void dequant_q5_0_block(const block_q5_0 *block, float *output)
Dequantize a single Q5_0 block to FP32.
Definition dequant_kernels.c:163
gemm_nt_q5_0_q8_0_m4n2
void gemm_nt_q5_0_q8_0_m4n2(const void *A_q8, const void *B_q5, const float *bias, float *C, int M, int N, int K)
Definition gemm_kernels_q5_0.c:2017
dequant_q5_0_row
void dequant_q5_0_row(const void *src, float *dst, size_t n_elements)
Dequantize Q5_0 row (multiple blocks)
Definition dequant_kernels.c:198
gemv_q5_0_ref
void gemv_q5_0_ref(float *y, const void *W, const float *x, int M, int K)
Matrix-vector multiply with Q5_0 weights (scalar reference)
Definition gemm_kernels_q5_0.c:67
gemv_q5_0_backward
void gemv_q5_0_backward(float *dX, const void *W, const float *dY, int M, int K)
Auto-dispatch backward.
Definition gemm_kernels_q5_0.c:753
gemv_q5_0_q8_0
void gemv_q5_0_q8_0(float *y, const void *W, const void *x_q8, int M, int K)
Matrix-vector multiply with Q5_0 weights and Q8_0 input.
Definition gemm_kernels_q5_0.c:1640
gemv_q5_0
void gemv_q5_0(float *y, const void *W, const float *x, int M, int K)
Auto-dispatch GEMV for Q5_0 weights based on CPU features.
Definition gemm_kernels_q5_0.c:550
vec_dot_q5_0_q8_0_ref
void vec_dot_q5_0_q8_0_ref(int n, float *s, const void *vx, const void *vy)
Quantized dot product: Q5_0 weights x Q8_0 input (scalar reference)
Definition gemm_kernels_q5_0.c:901
gemm_nt_q5_0
void gemm_nt_q5_0(const float *A, const void *B, const float *bias, float *C, int M, int N, int K)
Definition gemm_kernels_q5_0.c:833
vec_dot_q5_0_q8_0
void vec_dot_q5_0_q8_0(int n, float *s, const void *vx, const void *vy)
Auto-dispatch quantized dot product Q5_0 x Q8_0.
Definition gemm_kernels_q5_0.c:1602
gemm_q5_0
void gemm_q5_0(float *Y, const void *W, const float *X, int M, int N, int K)
Matrix-matrix multiply with Q5_0 weights.
Definition gemm_kernels_q5_0.c:684
gemm_nt_q5_0_q8_0_m4n2_tile
void gemm_nt_q5_0_q8_0_m4n2_tile(const void *A_q8, const void *B_q5, const float *bias, float *C, int M, int N, int K, int ldc)
Definition gemm_kernels_q5_0.c:1908
gemm_nt_q5_0_q8_0_m2n4_tile
void gemm_nt_q5_0_q8_0_m2n4_tile(const void *A_q8, const void *B_q5, const float *bias, float *C, int M, int N, int K, int ldc)
Definition gemm_kernels_q5_0.c:1766
gemm_nt_q5_0_ref
void gemm_nt_q5_0_ref(const float *A, const void *B, const float *bias, float *C, int M, int N, int K)
GEMM with transposed Q5_0 weights: C = A @ B^T.
Definition gemm_kernels_q5_0.c:790
gemm_nt_q5_0_q8_0
void gemm_nt_q5_0_q8_0(const void *A_q8, const void *B_q5, const float *bias, float *C, int M, int N, int K)
Batch GEMM with Q5_0 weights and Q8_0 activations for prefill.
Definition gemm_kernels_q5_0.c:1728
gemv_q5_0_parallel
void gemv_q5_0_parallel(float *y, const void *W, const float *x, int M, int K, int ith, int nth)
Parallel reference GEMV for Q5_0 × FP32.
Definition gemm_kernels_q5_0.c:579
dot_q5_0
float dot_q5_0(const void *w_q5_0, const float *x, int K)
Definition gemm_kernels_q5_0.c:872
gemm_nt_q5_0_q8_0_m2n4
void gemm_nt_q5_0_q8_0_m2n4(const void *A_q8, const void *B_q5, const float *bias, float *C, int M, int N, int K)
Definition gemm_kernels_q5_0.c:1894
gemv_q5_0_q8_0_parallel_simd
void gemv_q5_0_q8_0_parallel_simd(float *y, const void *W, const void *x_q8, int M, int K, int ith, int nth)
Parallel SIMD GEMV for Q5_0 x Q8_0 with prefetching.
Definition gemm_kernels_q5_0.c:1662
gemv_q5_0_parallel_simd
void gemv_q5_0_parallel_simd(float *y, const void *W, const float *x, int M, int K, int ith, int nth)
Parallel SIMD GEMV for Q5_0 × FP32 with prefetching.
Definition gemm_kernels_q5_0.c:625
C
#define C(color)
Definition show_config.c:39
block_q5_0
Definition ckernel_quant.h:70
block_q5_0::d
ck_half d
Definition ckernel_quant.h:71
block_q5_0::qh
uint8_t qh[4]
Definition ckernel_quant.h:72
block_q5_0::qs
uint8_t qs[32/2]
Definition ckernel_quant.h:73
block_q8_0
Definition ckernel_quant.h:104
block_q8_0::qs
int8_t qs[32]
Definition ckernel_quant.h:106
mask
int32_t int32_t int32_t int32_t int32_t mask
Definition tokenizer.h:234