Add Q4_1 tiled VNNI GEMM for prefill: pp256 +34%, pp2048 +42% vs vec_dot patch

README.md

@@ -36,10 +36,10 @@ Measured on Intel Xeon Platinum 8559C, 16 threads, with the patch (median of 3 i
 | Q4_K_XL, stock llama.cpp | 231 | 274 | 38.2 | 18.363 +- 0.795 |
 | **Q4_K_XL, patched** | **296** | **307** | 42.8 | 18.277 +- 0.787 |
 | Q4_1, stock llama.cpp | 183 | 189 | 48.3 | 18.260 +- 0.780 |
-| **Q4_1, patched** | 199 | 209 | **54.0** | 18.240 +- 0.779 |
+| **Q4_1, patched** | 266 | 297 | **55.7** | 18.240 +- 0.779 |
 
-- **Long prompts / large contexts (prefill-heavy): Q4_K_XL**. Its K-quant weights use the tiled VNNI matmul; patched prefill is ~47 % faster than Q4_1.
-- **Chat and agent loops with a cached prefix (decode-heavy): Q4_1**. It reads fewer bytes per weight (0.625) and decodes ~26 % faster.
+- **Long prompts / large contexts (prefill-heavy): Q4_K_XL** still prefills a little faster, but patched Q4_1 is now within ~10 % (it uses the tiled VNNI matmul too).
+- **Chat and agent loops with a cached prefix (decode-heavy): Q4_1**. It reads fewer bytes per weight (0.625) and decodes ~30 % faster.
 - Quality is equivalent within measurement error: the PPL differences above are smaller than the +-0.8 standard error of a 32-chunk test.
 - Single runs on this shared host vary by about +-8 %; the numbers are medians, raw data in `benchmarks/`.
 
@@ -48,7 +48,7 @@ Measured on Intel Xeon Platinum 8559C, 16 threads, with the patch (median of 3 i
 | file | pp256 | pp2048 | tg128 |
 |---|---|---|---|
 | Q4_K_XL | 1.28x | 1.12x | 1.12x |
-| Q4_1 | 1.09x | 1.10x | 1.12x |
+| Q4_1 | 1.46x | 1.57x | 1.15x |
 
 For Q4_K_XL most of this comes from the Qwen3.6 patch kernels (Q4_K/Q5_K/Q6_K VNNI, AMX tweaks, MoE tile grouping) also
 applying to flint; against the Qwen3.6 patch alone the flint-specific part is small (pp256 +12 %, decode +2 %, pp2048 equal).
@@ -56,22 +56,25 @@ The flint-specific changes:
 
 1. **IQ4_XS x Q8_K `vec_dot` in AVX512-VNNI** (x86 had AVX2 only). Unsigned LUT (+128) with `vpdpbusd`, bias removed through `bsums`, L2 prefetch.
 2. **Q4_1 x Q8_1 `vec_dot` in AVX-512VL VNNI**. Two blocks per iteration, packed into one ymm with memory-source `vinserti128`, low/high nibbles each through one `vpdpbusd`, fp16 headers converted together. A 512-bit version using cross-lane shuffles was slower than stock and dropped.
-3. Q5_K 5th-bit merge with `vptestmb` + masked add (bit-identical to the previous version, no speed change).
+3. **Q4_1 in the tiled VNNI GEMM** for prefill. Q4_1 has no super-block scale, so the tiles gain per-32-block f32 scale/min arrays (`sd`/`sm`) and a `Q41` kernel flag; activations reuse the existing Q8_1 `vec_dot_type`. Per 32-block: `d_a * (d_w * dot(q, code) + m_w * sum(code))`, integer dot/sum exact. This lifts pp256/pp2048 by ~34 %/42 % over the `vec_dot` path.
+4. Q5_K 5th-bit merge with `vptestmb` + masked add (bit-identical to the previous version, no speed change).
 
 ## Accuracy
 
-Wikitext-2 test, 32 chunks x 512, KL divergence against stock llama.cpp logits of the same file:
+Wikitext-2 test, 32 chunks x 512, KL divergence against the stated baseline logits of the same file:
 
-| file | KLD patched vs stock | same top-1 | noise floor: stock vs stock with `-fa off` |
+| comparison | KLD | same top-1 | noise floor: stock vs stock with `-fa off` |
 |---|---|---|---|
-| Q4_K_XL | 0.0666 +- 0.0049 | ~92 % | 0.0762 (also `-ub 128`: 0.0700) |
-| Q4_1 | 0.0364 +- 0.0029 | 93.8 % | 0.0385 |
+| Q4_K_XL, patched vs stock | 0.0666 +- 0.0049 | ~92 % | 0.0762 (also `-ub 128`: 0.0700) |
+| Q4_1, `vec_dot` patch vs stock | 0.0364 +- 0.0029 | 93.8 % | 0.0385 |
+| Q4_1, tiled GEMM vs `vec_dot` patch | 0.0372 +- 0.0026 | 94.2 % | 0.0385 |
 
 The KLD is **inside the noise floor of stock itself**: flint is very sensitive to tiny numeric differences (changing only the
 attention path of stock gives the same KLD), so it does not indicate a defect. Supporting checks: a graph-level `mul_mat` /
 `mul_mat_id` test gives the same relative L2 error versus a float reference as stock (Q4_K/Q5_K/Q6_K/IQ4_XS ~7e-3, Q4_1 ~5.4e-3, Q8_0
-~5.4e-3). Integer dot products are exact; only float accumulation order changes. Why small differences get amplified was not
-verified (top-6-of-256 expert routing ties is a guess). The earlier version of this repo reported KLD 0.2058 for a different patch.
+~5.4e-3; the tiled Q4_1 path gives the same 5.29e-3/5.35e-3 as the Q4_1 `vec_dot` path). Integer dot products are exact; only
+float accumulation order changes. Why small differences get amplified was not verified (top-6-of-256 expert routing ties is a
+guess). The earlier version of this repo reported KLD 0.2058 for a different patch.
 
 ## Where the time goes (why further kernel work is limited)
 
@@ -79,9 +82,10 @@ Measured DRAM read bandwidth on this host is 125-175 GB/s with 16 threads (singl
 weights per token (~1.6 GB dense Q8_0 attention/shared-expert/lm_head, ~0.5 GB for the 6 selected experts), i.e. ~90 GB/s effective:
 mostly memory-bound. Per-type effective bandwidth during decode (Q4_K_XL): Q8_0 dense ~100-170 GB/s (lm_head ~170), IQ4_XS experts ~100,
 Q5_K experts ~50. The Q5_K gap is not kernel throughput (cache-resident Q5_K and IQ4_XS run the same, 232 vs 218 GB/s); a
-prefetch-all-experts experiment did not help and was dropped; the cause was not found. Prefill is compute-bound; FlashAttention
-(24 %) and MoE matmul are the largest items. Q4_1 prefill still uses the per-row `vec_dot` path (60 % of prefill time is
-`MUL_MAT_ID q4_1`); adding it to the tiled GEMM would need per-32-block float scales in the tiled kernel and was not done.
+prefetch-all-experts experiment did not help and was dropped; the cause was not found. Prefill is compute-bound; once Q4_1 is
+in the tiled GEMM, FlashAttention and the tiled MoE matmul are the largest items (the per-row Q4_1 `vec_dot` remains only on the
+decode path). The Q4_1 tiled kernel carries per-32-block f32 scales/mins in the tiles (`sd`/`sm`) on top of the integer
+sub-block layout used by the K-quants.
 
 ## Build and run
 
@@ -96,13 +100,13 @@ values). `patches/optional-profiler.patch` adds `GGML_PROF=1` per-op timing.
 ## Limitations
 
 - Validated only on Intel Xeon Platinum 8559C (AVX512-VNNI + AMX), 16 threads, on a shared host. Other CPUs are untested; the new kernels need AVX512-VNNI (the Q4_1 one also AVX512VL/F16C).
-- The Q4_K_XL numbers were measured before the Q4_1 kernel was added; the K-quant code path is unchanged except for removing unused constants, and was not re-benchmarked.
+- The Q4_K_XL numbers were measured before the Q4_1 kernels were added; the K-quant code path is unchanged except for removing unused constants, and was not re-benchmarked (the Q4_1 tiled kernel was measured after).
 - PPL and KLD use a short test (32 x 512 tokens); no downstream task evaluation was run.
 
 ## Files
 
 - `patches/flint-cpu-opt.patch` (base `patches/BASE_COMMIT`), `patches/optional-profiler.patch`, `build.sh`, `machine.txt`
-- `benchmarks/summary.md` (Q4_K_XL), `benchmarks/q41-summary.md` (Q4_1), raw A/B CSVs, KLD logs, op profiles, KLD bisect, `graph_accuracy_test.cpp`
+- `benchmarks/summary.md` (Q4_K_XL), `benchmarks/q41-summary.md` (Q4_1), raw A/B CSVs (`ab-interleaved-flint.csv`, `q41-ab-interleaved.csv`, `q41-tiled-ab-interleaved.csv`), KLD logs (`kld_vs_stock.log`, `q41-kld.log`, `q41-tiled-kld.log`), op profiles, KLD bisect, `graph_accuracy_test.cpp`
 
 ## Credits and license
 

benchmarks/graph_accuracy_test.cpp

@@ -10,7 +10,7 @@
 int main(){
     ggml_backend_t be = ggml_backend_cpu_init(); ggml_backend_cpu_set_n_threads(be,16);
     std::mt19937 g(7); std::normal_distribution<float> nd(0,1);
-    ggml_type types[]={GGML_TYPE_Q4_K,GGML_TYPE_Q5_K,GGML_TYPE_Q6_K,GGML_TYPE_IQ4_XS,GGML_TYPE_Q8_0};
+    ggml_type types[]={GGML_TYPE_Q4_K,GGML_TYPE_Q5_K,GGML_TYPE_Q6_K,GGML_TYPE_IQ4_XS,GGML_TYPE_Q4_1,GGML_TYPE_Q8_0};
     const int K=2048,N=512,E=64,NU=6;
     for(ggml_type t:types) for(int mode=0;mode<2;mode++) for(int T: {1,64}){
         ggml_quantize_init(t);

benchmarks/q41-summary.md

@@ -2,18 +2,31 @@
 
 The Q4_1 file (24.2 GiB, every weight tensor Q4_1; token_embd and the router are F32) is not hosted here; download it from
 [RecursiveSI/flint-int4-gguf](https://huggingface.co/RecursiveSI/flint-int4-gguf). The patch adds an AVX512-VNNI `vec_dot` for
-Q4_1 x Q8_1 (the x86 path had AVX2 only).
+Q4_1 x Q8_1 (the x86 path had AVX2 only) and puts Q4_1 into the tiled VNNI GEMM used for prefill.
 
-`llama-bench -t 16 -fa on -lm mlock -b 1024 -ub 1024`, 3 interleaved rounds, median t/s (raw: `q41-ab-interleaved.csv`),
+`llama-bench -t 16 -fa on -lm mlock -b 1024 -ub 1024`, 3 interleaved rounds, median t/s (raw: `q41-ab-interleaved.csv` for the
+vec_dot kernel, `q41-tiled-ab-interleaved.csv` for the tiled GEMM),
 Xeon Platinum 8559C, same host as the other benchmarks (shared, +-8 % per run).
 
 | | pp256 | pp2048 | tg128 |
 |---|---|---|---|
 | stock llama.cpp (de5122d) | 182.7 | 188.7 | 48.3 |
-| this patch | **198.8** | **208.5** | **54.0** |
-| gain | +9 % | +10 % | +12 % |
-
-## Kernel
+| this patch, vec_dot only | 198.8 | 208.5 | 54.0 |
+| this patch, tiled GEMM | **266.5** | **297.2** | **55.7** |
+| gain vs stock | +46 % | +57 % | +15 % |
+| gain vs vec_dot only | +34 % | +42 % | +3 % (noise) |
+
+Decode is unchanged by the tiled path (single-row work stays below the tiled row threshold); the small tg difference is noise.
+
+## Tiled Q4_1 GEMM (prefill)
+Q4_1 tiles carry no super-block scale: each 32-wide block of a 256-element slab has its own f32 `d` (scale) and `m` (min).
+The tiled kernel gains `sd`/`sm` arrays (per-32-block f32, indexed per tile row) and a `Q41` template flag; activations use the
+existing Q8_1 `vec_dot_type`, so no new quantizer was needed. Per 32-block the kernel computes
+`d_a * (d_w * dot(q, code) + m_w * sum(code))`: the integer dot and integer code sum are exact and only the f32 accumulation
+order differs from `vec_dot`. VNNI uses `vpdpbusd` with the 0..15 weight nibbles as the unsigned side. The driver scales the
+k-block byte offset by `256/blck_size` (8 for Q4_1) instead of assuming one block per 256 elements.
+
+## vec_dot kernel (decode)
 Two 32-wide blocks per iteration. The two blocks' nibbles are packed into one ymm with memory-source `vinserti128`
 (no shuffle-port work), low and high nibbles each go through one `vpdpbusd` against the matching Q8_1 halves, the four fp16
 header values of the two blocks are converted at once, and the per-block `d*dy` scale is spread with a single `vpermps`.
@@ -21,22 +34,24 @@ A 512-bit version that re-ordered nibbles with cross-lane shuffles was slower th
 Prefetch distance `GGML_Q41_PF` (blocks, default 192; 192-320 equal in the cold-MoE microbenchmark).
 Microbenchmarks (16 threads, K=2048): cache-resident mul_mat 183 -> 220 GB/s, cold MoE decode (`mul_mat_id`) 43 -> 57 GB/s.
 
-## Accuracy (wikitext-2 test, 32 x 512, vs stock logits)
+## Accuracy
+wikitext-2 test, 32 x 512, tiled Q4_1 vs the vec_dot path of the same patch (logits baseline):
+
 | | KLD | same top-1 |
 |---|---|---|
-| this patch | 0.0364 +- 0.0029 | 93.8 % |
+| tiled Q4_1 vs vec_dot (`q41-tiled-kld.log`) | 0.0372 +- 0.0026 | 94.2 % |
+| vec_dot patch vs stock | 0.0364 +- 0.0029 | 93.8 % |
 | stock with `-fa off` (noise floor) | 0.0385 +- 0.0027 | 93.9 % |
 
-PPL 18.240 (patch) vs 18.260 (stock). Graph-level `mul_mat` / `mul_mat_id` error versus a float reference is identical to stock
-(5.3e-3 .. 5.4e-3). Integer dot products are exact; only float accumulation order differs.
+The tiled path is inside the patch's own noise floor. Graph-level `mul_mat` / `mul_mat_id` error versus a float reference is
+identical for tiled and vec_dot (`5.29e-3` vs `5.29e-3` for T=1, `5.35e-3` vs `5.39e-3` for T=64; Q8_0 reference 5.35e-3).
+Integer dot products are exact; only float accumulation order differs.
 
 ## Versus the Q4_K_XL file (same patch, same host)
 | file | pp256 | pp2048 | tg128 |
 |---|---|---|---|
 | flint-UD-Q4_K_XL (25.1 GB), patched | 296 | 307 | 42.8 |
-| flint-int4-q4_1 (26.0 GB), patched | 199 | 209 | 54.0 |
+| flint-int4-q4_1 (26.0 GB), patched | 266 | 297 | 55.7 |
 
-Q4_1 decodes about 26 % faster (0.625 B/weight and a simple layout) but prefills about 32 % slower: K-quant weights go through the
-tiled VNNI GEMM, while Q4_1 still uses the per-row `vec_dot` path (60 % of prefill time is `MUL_MAT_ID q4_1`). Supporting Q4_1 in
-the tiled path would need per-32-block float scales in the tiled kernel (its tiles carry integer sub-block scales under one
-super-block scale) and was not done. PPL of the two files is similar (18.26 vs 18.36 stock).
+Q4_1 decodes about 30 % faster (0.625 B/weight and a simple layout) and now prefills within about 10 % of Q4_K_XL (the tiled
+Q4_1 GEMM closed most of the earlier gap). PPL of the two files is similar (18.26 vs 18.36 stock).

benchmarks/q41-tiled-ab-interleaved.csv

@@ -0,0 +1,19 @@
+mode,round,test,t_s
+tiled_off,1,pp256,155.25
+tiled_off,1,pp2048,197.85
+tiled_off,1,tg128,55.74
+tiled_on,1,pp256,279.03
+tiled_on,1,pp2048,297.21
+tiled_on,1,tg128,51.19
+tiled_off,2,pp256,213.80
+tiled_off,2,pp2048,204.47
+tiled_off,2,tg128,55.76
+tiled_on,2,pp256,266.45
+tiled_on,2,pp2048,299.49
+tiled_on,2,tg128,55.67
+tiled_off,3,pp256,209.90
+tiled_off,3,pp2048,205.69
+tiled_off,3,tg128,58.09
+tiled_on,3,pp256,247.40
+tiled_on,3,pp2048,293.17
+tiled_on,3,tg128,56.45

benchmarks/q41-tiled-kld.log

@@ -0,0 +1,14 @@
+== base logits (this patch, Q4_1 prefill via vec_dot; GGML_CPU_TILED_MM=0)
+1.38.076.158 I Final estimate: PPL = 18.2396 +/- 0.77914
+
+== tiled Q4_1 vs vec_dot base (GGML_CPU_TILED_MM=1, --kl-divergence)
+chunk             PPL               ln(PPL(Q)/PPL(base))          KL Divergence              Δp RMS            Same top p
+Mean PPL(Q)                   :  18.272596 ±   0.781391
+Mean PPL(base)                :  17.736553 ±   0.734102
+Cor(ln(PPL(Q)), ln(PPL(base))):  98.98%
+Mean ln(PPL(Q)/PPL(base))     :   0.029775 ±   0.006166
+Mean    KLD:   0.037245 ±   0.002587
+Same top p: 94.154 ± 0.260 %
+
+wikitext-2 test, 32 chunks x 512, Xeon Platinum 8559C, 16 threads.
+For reference, the existing vec_dot patch is 0.036404 +- 0.002873 vs stock, and the stock -fa off noise floor is 0.038469 +- 0.002729.

patches/flint-cpu-opt.patch

@@ -757,11 +757,344 @@ index 9282ca299..0431b7ada 100644
              if (cur->ne[1] % 8 == 0) {
                  return &q4_K_8x8_q8_K;
              }
+diff --git a/ggml/src/ggml-cpu/tiled/tiled-kernel.cpp b/ggml/src/ggml-cpu/tiled/tiled-kernel.cpp
+index 5819ac07c..4d4eaae6f 100644
+--- a/ggml/src/ggml-cpu/tiled/tiled-kernel.cpp
++++ b/ggml/src/ggml-cpu/tiled/tiled-kernel.cpp
+@@ -204,6 +204,80 @@ static void tiled_run_microtile_vnni(const tiled_tile_src0 & src0, const tiled_t
+     tiled_run_micro_vnni_8x16<SUBBLK, HAS_MIN, BIAS, NK>(src0, src1, i0 + 8,  j0, num_k, slab, buf, buf_stride);
+ }
+ 
++// q4_1: every 32-wide block has its own f32 scale and min on both sides (weights: d,m; activations: d).
++// weight = d_w*q + m_w (q in 0..15), act = d_a*code (code int8), so one 32-block contributes
++//   d_a * (d_w * dot(q, code) + m_w * sum(code)).
++// The integer dot and the integer sum are exact; the four per-block f32 scale factors are applied in f32,
++// which only changes the accumulation order (same as the reference vec_dot, which also accumulates in f32).
++template <int NK>
++static void tiled_run_micro_q41_vnni_8x16(const tiled_tile_src0 & src0, const tiled_tile_src1 & src1,
++                                          int i0, int j0, int num_k, int slab, float * buf, int buf_stride) {
++    constexpr int SUBBLK = 32;
++    constexpr int NB = TILED_TILE_K / SUBBLK; // 8
++    constexpr int NS = SUBBLK / 16;           // 2
++    constexpr int NG = SUBBLK / 4;            // 8
++    constexpr int NUM_ROWS = 8;
++
++    const int nkr  = (NK > 0) ? NK : num_k;
++    const int seff = (NK > 0) ? 0 : slab;
++    const int qk_stride = nkr * TILED_TILE_K;
++    const int qk_off = seff * TILED_TILE_K;
++    const int nb_stride = NB * nkr;
++    const int nb_off = seff * NB;
++    const int bs_stride = (nkr == 1) ? TILED_TILE_ROWS : TILED_MICRO;
++
++    __m512 facc[NUM_ROWS];
++    for (int t = 0; t < NUM_ROWS; t++) { facc[t] = _mm512_setzero_ps(); }
++
++    const int32_t * bsums = src1.bsums;
++    for (int s = 0; s < NB; s++) {
++        __m512i bsums32 = _mm512_load_si512((const __m512i *) &bsums[(seff * TILED_MICRO + s * NS) * bs_stride + j0]);
++        for (int u = 1; u < NS; u++) {
++            bsums32 = _mm512_add_epi32(bsums32, _mm512_load_si512((const __m512i *) &bsums[(seff * TILED_MICRO + s * NS + u) * bs_stride + j0]));
++        }
++        const __m512 bterm = _mm512_cvtepi32_ps(bsums32);
++
++        __m512i acc16[NUM_ROWS];
++        for (int t = 0; t < NUM_ROWS; t++) { acc16[t] = _mm512_setzero_si512(); }
++
++#ifdef __GNUC__
++        #pragma GCC unroll 8
++#endif
++        for (int g = 0; g < NG; g++) {
++            const int kg = s * NG + g;
++            const __m512i codes = _mm512_load_si512((const __m512i *) &src1.q[(j0 / TILED_MICRO) * (TILED_MICRO * qk_stride)
++                + qk_off + (kg % TILED_MICRO) * qk_stride + (kg / TILED_MICRO) * (TILED_MICRO * 4)]);
++            for (int t = 0; t < NUM_ROWS; t++) {
++                const uint32_t u4 = *(const uint32_t *) &src0.q[(i0 + t) * qk_stride + qk_off + kg * 4];
++                const __m512i u4b = _mm512_set1_epi32((int) u4);
++                acc16[t] = _mm512_dpbusd_epi32(acc16[t], u4b, codes);
++            }
++        }
++
++        const __m512 dav = _mm512_loadu_ps(&src1.sd[(seff * NB + s) * bs_stride + j0]);
++        for (int t = 0; t < NUM_ROWS; t++) {
++            const int ar = i0 + t;
++            const float dw = src0.sd[ar * nb_stride + nb_off + s];
++            const float mw = src0.sm[ar * nb_stride + nb_off + s];
++            const __m512 raw = _mm512_cvtepi32_ps(acc16[t]);
++            const __m512 tv  = _mm512_fmadd_ps(_mm512_set1_ps(dw), raw, _mm512_mul_ps(_mm512_set1_ps(mw), bterm));
++            facc[t] = _mm512_fmadd_ps(tv, dav, facc[t]);
++        }
++    }
++
++    for (int t = 0; t < NUM_ROWS; t++) {
++        float * p = &buf[(i0 + t) * buf_stride + j0];
++        _mm512_store_ps(p, _mm512_add_ps(_mm512_load_ps(p), facc[t]));
++    }
++}
++
++template <int NK>
++static void tiled_run_microtile_q41_vnni(const tiled_tile_src0 & src0, const tiled_tile_src1 & src1,
++                                         int i0, int j0, int num_k, int slab, float * buf, int buf_stride) {
++    tiled_run_micro_q41_vnni_8x16<NK>(src0, src1, i0,      j0, num_k, slab, buf, buf_stride);
++    tiled_run_micro_q41_vnni_8x16<NK>(src0, src1, i0 + 8,  j0, num_k, slab, buf, buf_stride);
++}
++
+ #endif // __AVX512VNNI__ && __AVX512VL__
+ 
+ #if defined(__AVX2__)
+@@ -496,23 +570,79 @@ static void tiled_run_microtile_avx(const tiled_tile_src0 & src0, const tiled_ti
+ 
+ #endif // __AVX__ && !__AVX2__
+ 
++// q4_1 reference/scalar path (non-VNNI builds): same math as the VNNI kernel, per 32-block f32 scales.
++template <int NK>
++static void tiled_run_microtile_q41_scalar(const tiled_tile_src0 & src0, const tiled_tile_src1 & src1,
++                                           int i0, int j0, int num_k, int slab, float * buf, int buf_stride) {
++    constexpr int SUBBLK = 32;
++    constexpr int NB = TILED_TILE_K / SUBBLK;
++    constexpr int NS = SUBBLK / 16;
++
++    const int nkr  = (NK > 0) ? NK : num_k;
++    const int seff = (NK > 0) ? 0 : slab;
++    const int qk_stride = nkr * TILED_TILE_K;
++    const int qk_off = seff * TILED_TILE_K;
++    const int nb_stride = NB * nkr;
++    const int nb_off = seff * NB;
++    const int bs_stride = (nkr == 1) ? TILED_TILE_ROWS : TILED_MICRO;
++
++    for (int t = 0; t < TILED_MICRO; t++) {
++        const int ar = i0 + t;
++        const uint8_t * q0 = &src0.q[ar * qk_stride + qk_off];
++        for (int j = 0; j < TILED_MICRO; j++) {
++            const int br = j0 + j;
++            const int8_t * q1 = &src1.q[br * qk_stride + qk_off];
++            float acc = 0.0f;
++            for (int s = 0; s < NB; s++) {
++                int32_t bsum = 0;
++                for (int u = 0; u < NS; u++) { bsum += src1.bsums[(seff * TILED_MICRO + s * NS + u) * bs_stride + br]; }
++                int32_t raw = 0;
++                for (int e = 0; e < SUBBLK; e++) { raw += (int32_t) q0[s * SUBBLK + e] * (int32_t) q1[s * SUBBLK + e]; }
++                const float da = src1.sd[(seff * NB + s) * bs_stride + br];
++                const float dw = src0.sd[ar * nb_stride + nb_off + s];
++                const float mw = src0.sm[ar * nb_stride + nb_off + s];
++                acc += da * (dw * (float) raw + mw * (float) bsum);
++            }
++            buf[ar * buf_stride + j0 + j] += acc;
++        }
++    }
++}
++
+ // main microtile entry point. num_k = K-blocks per row (the narrow path reads num_k slabs at row
+ // stride num_k*256 so each weight row is one long stream); slab = this call's slab. num_k=1,
+ // slab=0 is the standard single-slab path.
+ // num_k==1 dispatches to the NK=1 instantiation, whose strides/offsets are compile-time
+ // constants (shifts/scaled-leas). num_k>1 (the narrow path) uses the NK=0 runtime-strided
+ // version; it is memory-bound, so one version serves all of it.
+-template <int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS>
++// Q41: q4_1 weights + q8_1 activations (per-32-block f32 scales) use the dedicated kernels above.
++template <int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS, bool Q41>
+ void tiled_run_microtile(const tiled_tile_src0 & src0, const tiled_tile_src1 & src1,
+                          int i0, int j0, int num_k, int slab, float * buf, int buf_stride) {
+     const bool standard = (num_k == 1);
+ #if defined(__AVX512VNNI__) && defined(__AVX512VL__) && defined(__AVX512DQ__)
++    if constexpr (Q41) {
++        if (standard) {
++            tiled_run_microtile_q41_vnni<1>(src0, src1, i0, j0, num_k, slab, buf, buf_stride);
++        } else {
++            tiled_run_microtile_q41_vnni<0>(src0, src1, i0, j0, num_k, slab, buf, buf_stride);
++        }
++        return;
++    }
+     if (standard) {
+         tiled_run_microtile_vnni<SUBBLK, HAS_MIN, BIAS, 1>(src0, src1, i0, j0, num_k, slab, buf, buf_stride);
+     } else {
+         tiled_run_microtile_vnni<SUBBLK, HAS_MIN, BIAS, 0>(src0, src1, i0, j0, num_k, slab, buf, buf_stride);
+     }
+-#elif defined(__AVX2__)
++#else
++    if constexpr (Q41) {
++        if (standard) {
++            tiled_run_microtile_q41_scalar<1>(src0, src1, i0, j0, num_k, slab, buf, buf_stride);
++        } else {
++            tiled_run_microtile_q41_scalar<0>(src0, src1, i0, j0, num_k, slab, buf, buf_stride);
++        }
++        return;
++    }
++#if defined(__AVX2__)
+     if (standard) {
+         tiled_run_microtile_avx2<SUBBLK, HAS_MIN, BIAS, 1, ACTBIAS>(src0, src1, i0, j0, num_k, slab, buf, buf_stride);
+     } else {
+@@ -531,6 +661,7 @@ void tiled_run_microtile(const tiled_tile_src0 & src0, const tiled_tile_src1 & s
+         tiled_run_microtile_scalar<SUBBLK, HAS_MIN, BIAS, 0>(src0, src1, i0, j0, num_k, slab, buf, buf_stride);
+     }
+ #endif
++#endif
+ }
+ 
+ // explicit instantiations for the in-use formats (q4_K and q5_K share the constants)
+@@ -556,6 +687,9 @@ template void tiled_run_microtile<16, false, 4, true>(const tiled_tile_src0 & sr
+ // q2_K: BIAS = 0
+ template void tiled_run_microtile<16, true, 0, false>(const tiled_tile_src0 & src0, const tiled_tile_src1 & src1,
+                                                int i0, int j0, int num_k, int slab, float * buf, int buf_stride);
++// q4_1: per-32-block f32 scales, q8_1 activations
++template void tiled_run_microtile<32, false, 0, false, true>(const tiled_tile_src0 & src0, const tiled_tile_src1 & src1,
++                                                   int i0, int j0, int num_k, int slab, float * buf, int buf_stride);
+ 
+ 
+ #define MIN(a, b) ((a) < (b) ? (a) : (b))
+diff --git a/ggml/src/ggml-cpu/tiled/tiled-kernel.h b/ggml/src/ggml-cpu/tiled/tiled-kernel.h
+index 1a97f406a..4bd9642cc 100644
+--- a/ggml/src/ggml-cpu/tiled/tiled-kernel.h
++++ b/ggml/src/ggml-cpu/tiled/tiled-kernel.h
+@@ -29,12 +29,16 @@
+ // SUBBLK=32 formats index at stride 8 and leave the slack unused.
+ struct tiled_tile_src0 {
+     static constexpr int NB_MAX = TILED_TILE_K / 16; // max subblocks per 256-elem block
++    static constexpr int NB32  = TILED_TILE_K / 32;  // 32-wide subblocks (q4_1)
+ 
+     alignas(64) uint8_t q[TILED_TILE_ROWS * TILED_TILE_K]; // unsigned quants, widened to uint8
+     float    d[TILED_TILE_ROWS];  // One d from each input block, widened to f32
+     float    dmin[TILED_TILE_ROWS]; // dmin from each input block (if applicable), widened to F32
+     int32_t   scales[TILED_TILE_ROWS * NB_MAX];  // per-subblock scale, stored as int32_t
+     int32_t   mins[TILED_TILE_ROWS * NB_MAX];    // per-subblock min, used when HAS_MIN
++    // q4_1: no super-block scale, each 32-wide block carries its own f32 scale/min
++    float    sd[TILED_TILE_ROWS * NB32];  // per-32-block scale
++    float    sm[TILED_TILE_ROWS * NB32];  // per-32-block min
+ };
+ 
+ // src1 tile: built from q8_K (wdata)
+@@ -45,6 +49,8 @@ struct tiled_tile_src1 {
+     alignas(64) int32_t bsums[(TILED_TILE_K / 16) * TILED_TILE_ROWS];
+     // f32 (not f16): q8_k stores fp16, the unpack converts once
+     float       d[TILED_TILE_ROWS];
++    // q4_1: q8_1 activations carry a scale per 32-wide block (no super-block d)
++    float       sd[(TILED_TILE_K / 32) * TILED_TILE_ROWS];
+ };
+ 
+ // per-thread workspace: all tiled state lives here, allocated in wdata (one slot per thread)
+@@ -71,6 +77,12 @@ inline void tiled_unpk_nib4(const uint8_t * src, uint8_t * lo, uint8_t * hi) {
+     _mm256_storeu_si256((__m256i *) lo, _mm256_and_si256(v, _mm256_set1_epi8(0x0F)));
+     _mm256_storeu_si256((__m256i *) hi, _mm256_srli_epi32(_mm256_and_si256(v, _mm256_set1_epi8((int8_t) 0xF0)), 4));
+ }
++// 16-byte source (one q4_1 block: lo -> elems 0..15, hi -> elems 16..31)
++inline void tiled_unpk_nib4_16(const uint8_t * src, uint8_t * lo, uint8_t * hi) {
++    const __m128i v = _mm_loadu_si128((const __m128i *) src);
++    _mm_storeu_si128((__m128i *) lo, _mm_and_si128(v, _mm_set1_epi8(0x0F)));
++    _mm_storeu_si128((__m128i *) hi, _mm_srli_epi32(_mm_and_si128(v, _mm_set1_epi8((int8_t) 0xF0)), 4));
++}
+ // 2-bit values at bit offset S
+ template <int S> inline void tiled_unpk_2bit(const uint8_t * src, uint8_t * dst) {
+     _mm256_storeu_si256((__m256i *) dst, _mm256_and_si256(
+@@ -129,6 +141,9 @@ inline void tiled_unpk_tern8(const uint8_t * src, int8_t delta, uint8_t * dst) {
+ inline void tiled_unpk_nib4(const uint8_t * src, uint8_t * lo, uint8_t * hi) {
+     for (int l = 0; l < 32; l++) { lo[l] = (uint8_t) (src[l] & 0xF); hi[l] = (uint8_t) (src[l] >> 4); }
+ }
++inline void tiled_unpk_nib4_16(const uint8_t * src, uint8_t * lo, uint8_t * hi) {
++    for (int l = 0; l < 16; l++) { lo[l] = (uint8_t) (src[l] & 0xF); hi[l] = (uint8_t) (src[l] >> 4); }
++}
+ template <int S>
+ inline void tiled_unpk_2bit(const uint8_t * src, uint8_t * dst) {
+     for (int l = 0; l < 32; l++) { dst[l] = (uint8_t) ((src[l] >> S) & 3); }
+@@ -164,7 +179,9 @@ inline void tiled_unpk_tern8(const uint8_t * src, int8_t delta, uint8_t * dst) {
+ // ACTBIAS (AVX2 only): the activation is pre-biased +128 by tiled_repack_src1,
+ // num_k = K-blocks per row: the tile holds num_k slabs at row stride num_k*256
+ // (Default case is num_k=1, 256x256 tiles, we go to longer num_k to improve memory bandwidth when num_rows is small)
+-template <int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS>
++// Q41: src0 is q4_1 (per-32-block f32 scale/min in sd/sm) and src1 is q8_1 (per-32-block f32 scale in sd);
++// the other format constants (SUBBLK=32, HAS_MIN=false, BIAS=0, ACTBIAS=false) are unused in that mode.
++template <int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS, bool Q41 = false>
+ void tiled_run_microtile(const tiled_tile_src0 & src0, const tiled_tile_src1 & src1,
+                          int i0, int j0, int num_k, int slab, float * buf, int buf_stride);
+ 
 diff --git a/ggml/src/ggml-cpu/tiled/tiled.cpp b/ggml/src/ggml-cpu/tiled/tiled.cpp
-index 77e1a9945..06fd9a2d5 100644
+index 77e1a9945..225ec0e78 100644
 --- a/ggml/src/ggml-cpu/tiled/tiled.cpp
 +++ b/ggml/src/ggml-cpu/tiled/tiled.cpp
-@@ -716,7 +716,9 @@ static void tiled_store_window(const float * buf, int n_src0, int n_src1, int bu
+@@ -534,11 +534,33 @@ static void tiled_unpack_src0(const block_iq1_m * rows, int64_t row_stride, int
+     }
+ }
+ 
++// q4_1: 8 x 32-wide blocks per 256-slab, each with its own f32 d (scale) and m (min).
++static void tiled_unpack_src0(const block_q4_1 * rows, int64_t row_stride, int n_rows, tiled_tile_src0 * tile, int num_k) {
++    GGML_ASSERT(n_rows <= TILED_TILE_ROWS);
++    constexpr int NB = TILED_TILE_K / 32; // 8 x 32-wide blocks
++
++    const int qk_stride = num_k * TILED_TILE_K;
++    const int nb_stride = NB * num_k;
++    for (int slab = 0; slab < num_k; slab++) {
++        for (int r = 0; r < n_rows; r++) {
++            const int q_off = r * qk_stride + slab * TILED_TILE_K;
++            const int s_off = r * nb_stride + slab * NB;
++            const block_q4_1 * x = rows + r * row_stride + slab * NB;
++            for (int s = 0; s < NB; s++) {
++                tile->sd[s_off + s] = ggml_fp16_to_fp32(x[s].data.data.d);
++                tile->sm[s_off + s] = ggml_fp16_to_fp32(x[s].data.data.m);
++                // 16 bytes -> elems 0..15 (low nibble) and 16..31 (high nibble)
++                tiled_unpk_nib4_16(x[s].qs, &tile->q[q_off + s * 32], &tile->q[q_off + s * 32 + 16]);
++            }
++        }
++    }
++}
++
+ // unpack src1 tile from q8_K rows. num_k = K-blocks per row: each row decodes num_k consecutive
+ // source blocks (rows[r][kblk + slab]) into the tile at row stride num_k*256, so each activation
+ // row is one long stream. num_k=1 is the standard single-slab path (byte-identical to before).
+-static void tiled_unpack_src1_q8_K(const block_q8_K * const * rows, int n_rows, tiled_tile_src1 * tile,
+-                                   int kblk, int num_k) {
++static void tiled_unpack_src1(const block_q8_K * const * rows, int n_rows, tiled_tile_src1 * tile,
++                              int kblk, int num_k) {
+     GGML_ASSERT(n_rows <= TILED_TILE_ROWS);
+     const int n_padded = (n_rows + TILED_MICRO - 1) & ~(TILED_MICRO - 1);
+     const int bs_stride = (num_k == 1) ? TILED_TILE_ROWS : TILED_MICRO;
+@@ -569,6 +591,39 @@ static void tiled_unpack_src1_q8_K(const block_q8_K * const * rows, int n_rows,
+     }
+ }
+ 
++// q4_1 src1: q8_1 activation rows. 8 blocks of 32 per 256-slab, each block has its own f32 scale;
++// the integer per-16 sums are recomputed from the codes (q8_1 only stores d and d*sum).
++static void tiled_unpack_src1(const block_q8_1 * const * rows, int n_rows, tiled_tile_src1 * tile,
++                              int kblk, int num_k) {
++    GGML_ASSERT(n_rows <= TILED_TILE_ROWS);
++    const int n_padded = (n_rows + TILED_MICRO - 1) & ~(TILED_MICRO - 1);
++    const int bs_stride = (num_k == 1) ? TILED_TILE_ROWS : TILED_MICRO;
++    constexpr int NBLK = TILED_TILE_K / 32; // 8 q8_1 blocks per 256-slab
++    constexpr int n16 = TILED_TILE_K / TILED_MICRO;
++    for (int slab = 0; slab < num_k; slab++) {
++        const int q_off = slab * TILED_TILE_K;
++        for (int r = 0; r < n_padded; r++) {
++            if (r < n_rows) {
++                const block_q8_1 * xb = rows[r] + (kblk + slab) * NBLK;
++                for (int b = 0; b < NBLK; b++) {
++                    memcpy(&tile->q[r * (num_k * TILED_TILE_K) + q_off + b * 32], xb[b].qs, 32);
++                    int32_t s0 = 0, s1 = 0;
++                    for (int e = 0; e < 16; e++) { s0 += xb[b].qs[e]; s1 += xb[b].qs[16 + e]; }
++                    tile->bsums[(slab * n16 + 2 * b + 0) * bs_stride + r] = s0;
++                    tile->bsums[(slab * n16 + 2 * b + 1) * bs_stride + r] = s1;
++                    tile->sd[(slab * NBLK + b) * bs_stride + r] = ggml_fp16_to_fp32(xb[b].data.data.d);
++                }
++                tile->d[slab * TILED_MICRO + r] = 1.0f; // per-32 scales are applied in the kernel
++            } else {
++                memset(&tile->q[r * (num_k * TILED_TILE_K) + q_off], 0, TILED_TILE_K);
++                for (int i = 0; i < n16; i++) { tile->bsums[(slab * n16 + i) * bs_stride + r] = 0; }
++                for (int b = 0; b < NBLK; b++) { tile->sd[(slab * NBLK + b) * bs_stride + r] = 0.0f; }
++                tile->d[slab * TILED_MICRO + r] = 1.0f;
++            }
++        }
++    }
++}
++
+ // GGML_CPU_TILED_MM: master switch, on by default. If off, we fast return false and normal vec_dot mul_mat resumes
+ static bool ggml_tiled_matmul_enabled(void) {
+     static bool enabled = true;
+@@ -617,6 +672,7 @@ static bool ggml_tiled_supported(const struct ggml_tensor * src0,
+         case GGML_TYPE_Q6_K:
+         case GGML_TYPE_Q5_K:
+         case GGML_TYPE_Q4_K:
++        case GGML_TYPE_Q4_1:
+         case GGML_TYPE_Q3_K:
+         case GGML_TYPE_Q2_K:
+         case GGML_TYPE_IQ4_XS:
+@@ -716,7 +772,9 @@ static void tiled_store_window(const float * buf, int n_src0, int n_src1, int bu
  // unpack and microtile are reused unchanged from mul_mat.
  
  // src0 rows per g group; finer than the TILED_TILE_K kernel tile so all threads stay busy at small ne01
@@ -772,6 +1105,218 @@ index 77e1a9945..06fd9a2d5 100644
  
  // like tiled_store_window, but the dst columns are not contiguous: col_ptrs[j] points at the
  // start of dst column j, whose rows are contiguous (dim 0, nb0 == 4 bytes).
+@@ -756,11 +814,11 @@ static void tiled_store_window_scatter(const float * buf, int n_src0, int n_src1
+ 
+ // one (g, k) macrotile: zero the acc window, sweep K in 256 element slabs, scatter the result rows
+ // rows points at this k window's routed src1 rows, row r at its base block
+-template <typename B, int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS>
++template <typename B, int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS, typename BA = block_q8_K, bool Q41 = false>
+ static void tiled_mmid_gemm_window(struct ggml_tensor * dst, const struct ggml_tensor * src0,
+                                    const char * src0_cur, int64_t r, int64_t k, int64_t nrows,
+                                    const int32_t * expert_rows,
+-                                   const block_q8_K * const * rows,
++                                   const BA * const * rows,
+                                    tiled_ws * ws) {
+     const int64_t ne00 = src0->ne[0];
+     const int64_t ne01 = src0->ne[1];
+@@ -770,6 +828,8 @@ static void tiled_mmid_gemm_window(struct ggml_tensor * dst, const struct ggml_t
+ 
+     const size_t src0_bs = ggml_type_size(src0->type);
+     const int64_t src0_stride = src0->nb[1] / src0_bs;
++    // source blocks per 256-element slab (1 for K-quants, 8 for q4_1's 32-wide blocks)
++    const int kblk_scale = (int) (TILED_TILE_K / ggml_blck_size(src0->type));
+ 
+     // the window is at most TILED_MMID_GROUP x TILED_TILE_K, so zero only that region of acc
+     for (int64_t i = 0; i < n_src0; i++) {
+@@ -792,19 +852,19 @@ static void tiled_mmid_gemm_window(struct ggml_tensor * dst, const struct ggml_t
+             const int kstart = (int) (k0 / TILED_TILE_K);
+             // activation: one 16-row band (nrows <= 16). Unpack the whole k_extent chunk (all
+             // slabs) so each activation row is a long stream, then repack each slab in place
+-            tiled_unpack_src1_q8_K(rows, nrows, &ws->src1, kstart, num_k);
++            tiled_unpack_src1(rows, nrows, &ws->src1, kstart, num_k);
+             // repack the whole k_extent chunk in place: num_k slabs x 4 tiles, 16 rows, row stride num_k*256
+             tiled_repack_src1(&ws->src1, 0, num_k, ACTBIAS);
+             // weight groups (16 at a time): unpack the k_extent chunk of the group (the long per-row
+             // read), then one standard MAC per slab accumulating into the same acc rows
+             for (int64_t ir0 = r; ir0 < r_end; ir0 += TILED_MICRO) {
+                 const int n0 = (int) MIN(TILED_MICRO, r_end - ir0);
+-                const B * wbase = (const B *) (src0_cur + ir0 * src0->nb[1] + kstart * src0_bs);
++                const B * wbase = (const B *) (src0_cur + ir0 * src0->nb[1] + kstart * kblk_scale * src0_bs);
+                 tiled_unpack_src0(wbase, src0_stride, n0, &ws->src0, num_k);
+                 tiled_repack_src0<SUBBLK>(&ws->src0, n0, num_k, BIAS, ACTBIAS);
+                 float * buf = ws->acc + (ir0 - r) * TILED_TILE_ROWS;
+                 for (int slab = 0; slab < num_k; slab++) {
+-                    tiled_run_microtile<SUBBLK, HAS_MIN, BIAS, ACTBIAS>(ws->src0, ws->src1, 0, 0, num_k, slab, buf, TILED_TILE_ROWS);
++                    tiled_run_microtile<SUBBLK, HAS_MIN, BIAS, ACTBIAS, Q41>(ws->src0, ws->src1, 0, 0, num_k, slab, buf, TILED_TILE_ROWS);
+                 }
+             }
+         }
+@@ -812,15 +872,15 @@ static void tiled_mmid_gemm_window(struct ggml_tensor * dst, const struct ggml_t
+         // standard per-slab path: K is stepped in 256-element slabs
+         for (int64_t ib = 0; ib < ne00; ib += TILED_TILE_K) {
+             const int kblk = (int) (ib / TILED_TILE_K);
+-            tiled_unpack_src0((const B *) (src0_cur + r * src0->nb[1] + kblk * src0_bs), src0_stride, n_src0, &ws->src0, 1);
++            tiled_unpack_src0((const B *) (src0_cur + r * src0->nb[1] + kblk * kblk_scale * src0_bs), src0_stride, n_src0, &ws->src0, 1);
+             tiled_repack_src0<SUBBLK>(&ws->src0, n_src0, 1, BIAS, ACTBIAS);
+-            tiled_unpack_src1_q8_K(rows, nrows, &ws->src1, kblk, 1);
++            tiled_unpack_src1(rows, nrows, &ws->src1, kblk, 1);
+             // 16x16 microtiles sweeping the window; repack each src1 band just-in-time
+             // (j0-outer) so only the bands actually used are repacked and each is L1-hot
+             for (int64_t ir1 = 0; ir1 < nrows; ir1 += TILED_MICRO) {
+                 tiled_repack_src1(&ws->src1, (int) ir1, 1, ACTBIAS);
+                 for (int64_t ir0 = r; ir0 < r_end; ir0 += TILED_MICRO) {
+-                    tiled_run_microtile<SUBBLK, HAS_MIN, BIAS, ACTBIAS>(ws->src0, ws->src1,
++                    tiled_run_microtile<SUBBLK, HAS_MIN, BIAS, ACTBIAS, Q41>(ws->src0, ws->src1,
+                         (int) (ir0 - r), (int) ir1, 1, 0,
+                         ws->acc, TILED_TILE_ROWS);
+                 }
+@@ -834,7 +894,7 @@ static void tiled_mmid_gemm_window(struct ggml_tensor * dst, const struct ggml_t
+ // one expert of MUL_MAT_ID: each k window's dispatched rows are pointed at by a per-window
+ // row pointer list and swept over the thread's row windows, the dst rows are scattered back
+ // per window
+-template <typename B, int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS>
++template <typename B, int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS, typename BA = block_q8_K, bool Q41 = false>
+ static void ggml_compute_forward_mul_mat_id_tiled_one_expert(
+         const struct ggml_compute_params * params,
+               struct ggml_tensor *         dst,
+@@ -877,22 +937,22 @@ static void ggml_compute_forward_mul_mat_id_tiled_one_expert(
+     for (int64_t k = 0; k < cne1; k += TILED_TILE_K) {
+         const int64_t nrows = MIN(TILED_TILE_K, cne1 - k);
+ 
+-        const block_q8_K * rows[TILED_TILE_ROWS];
++        const BA * rows[TILED_TILE_ROWS];
+         for (int64_t i = 0; i < nrows; i++) {
+             const int64_t i11 = expert_rows[2 * (k + i) + 0] % ne11;
+             const int64_t i12 = expert_rows[2 * (k + i) + 1];
+-            rows[i] = (const block_q8_K *) ((const char *) wdata + (i11 + i12 * ne11) * nbw1);
++            rows[i] = (const BA *) ((const char *) wdata + (i11 + i12 * ne11) * nbw1);
+         }
+ 
+         for (int64_t g = g0; g < g1; g++) {
+             const int64_t r = g * TILED_MMID_GROUP;
+ 
+-            tiled_mmid_gemm_window<B, SUBBLK, HAS_MIN, BIAS, ACTBIAS>(dst, src0, src0_cur, r, k, nrows, expert_rows, rows, ws);
++            tiled_mmid_gemm_window<B, SUBBLK, HAS_MIN, BIAS, ACTBIAS, BA, Q41>(dst, src0, src0_cur, r, k, nrows, expert_rows, rows, ws);
+         }
+     }
+ }
+ 
+-template <typename B, int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS>
++template <typename B, int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS, typename BA = block_q8_K, bool Q41 = false>
+ static void ggml_compute_forward_mul_mat_tiled_one_chunk(
+     const struct ggml_compute_params * params,
+     struct ggml_tensor * dst,
+@@ -928,6 +988,8 @@ static void ggml_compute_forward_mul_mat_tiled_one_chunk(
+ 
+     const int64_t src0_stride = nb01 / src0_bs;  // blocks between src0 rows
+     const int64_t src1_stride = (src1->type == vec_dot_type ? src1->nb[1] : row_size) / src1_bs;
++    // source blocks per 256-element slab (1 for K-quants, 8 for q4_1's 32-wide blocks)
++    const int kblk_scale = (int) (TILED_TILE_K / ggml_blck_size(src0->type));
+ 
+     const int64_t TILE = 256;
+     const int64_t MICRO = 16;
+@@ -960,9 +1022,9 @@ static void ggml_compute_forward_mul_mat_tiled_one_chunk(
+         char * dst_col = (char *) dst->data + i12 * nb2 + i13 * nb3;
+ 
+         // rows[r] is the window's row r base block; the k-slab offset is applied in the unpack
+-        const block_q8_K * rows[TILED_TILE_ROWS];
++        const BA * rows[TILED_TILE_ROWS];
+         for (int r = 0; r < n_src1; r++) {
+-            rows[r] = (const block_q8_K *) ((const char *) wdata + (iir1 + r) * src1_stride * src1_bs);
++            rows[r] = (const BA *) ((const char *) wdata + (iir1 + r) * src1_stride * src1_bs);
+         }
+ 
+         for (int64_t iir0 = ir0_start; iir0 < ir0_end; iir0 += TILE) {
+@@ -981,19 +1043,19 @@ static void ggml_compute_forward_mul_mat_tiled_one_chunk(
+                     const int kstart = (int) (k0 / TILED_TILE_K);
+                     // activation: one 16-row band (n_src1 <= 16). Unpack the whole k_extent chunk
+                     // (all slabs) so each activation row is a long stream, then repack each slab in place
+-                    tiled_unpack_src1_q8_K(rows, n_src1, &ws->src1, kstart, num_k);
++                    tiled_unpack_src1(rows, n_src1, &ws->src1, kstart, num_k);
+                     // repack the whole k_extent chunk in place, if kernel wants to: num_k slabs x 4 tiles, 16 rows, row stride num_k*256
+                     tiled_repack_src1(&ws->src1, 0, num_k, ACTBIAS);
+                     // weight groups: unpack the k_extent chunk of each 16-row group (the long per-row
+                     // read), then one standard MAC per slab accumulating into the same acc rows
+                     for (int64_t ir0 = iir0; ir0 < iir0_end; ir0 += MICRO) {
+                         const int n0 = (int) MIN(MICRO, iir0_end - ir0);
+-                        const B * wbase = (const B *) (src0_row + ir0 * nb01 + kstart * src0_bs);
++                        const B * wbase = (const B *) (src0_row + ir0 * nb01 + kstart * kblk_scale * src0_bs);
+                         tiled_unpack_src0(wbase, src0_stride, n0, &ws->src0, num_k);
+                         tiled_repack_src0<SUBBLK>(&ws->src0, n0, num_k, BIAS, ACTBIAS);
+                         float * buf = ws->acc + (ir0 - iir0) * TILED_TILE_ROWS;
+                         for (int slab = 0; slab < num_k; slab++) {
+-                            tiled_run_microtile<SUBBLK, HAS_MIN, BIAS, ACTBIAS>(ws->src0, ws->src1, 0, 0, num_k, slab, buf, TILED_TILE_ROWS);
++                            tiled_run_microtile<SUBBLK, HAS_MIN, BIAS, ACTBIAS, Q41>(ws->src0, ws->src1, 0, 0, num_k, slab, buf, TILED_TILE_ROWS);
+                         }
+                     }
+                 }
+@@ -1001,15 +1063,15 @@ static void ggml_compute_forward_mul_mat_tiled_one_chunk(
+                 // standard per-slab path: K is stepped in 256-element chunks
+                 for (int64_t ib = 0; ib < ne00; ib += TILE) {
+                     const int kblk = (int) (ib / TILE);
+-                    tiled_unpack_src0((const B *) (src0_row + iir0 * nb01 + kblk * src0_bs), src0_stride, n_src0, &ws->src0, 1);
++                    tiled_unpack_src0((const B *) (src0_row + iir0 * nb01 + kblk * kblk_scale * src0_bs), src0_stride, n_src0, &ws->src0, 1);
+                     tiled_repack_src0<SUBBLK>(&ws->src0, n_src0, 1, BIAS, ACTBIAS);
+-                    tiled_unpack_src1_q8_K(rows, n_src1, &ws->src1, kblk, 1);
++                    tiled_unpack_src1(rows, n_src1, &ws->src1, kblk, 1);
+ 
+                     // 16x16 microtiles sweeping the window; repack each src1 band before first use
+                     for (int64_t ir1 = iir1; ir1 < iir1_end; ir1 += MICRO) {
+                         tiled_repack_src1(&ws->src1, (int) (ir1 - iir1), 1, ACTBIAS);
+                         for (int64_t ir0 = iir0; ir0 < iir0_end; ir0 += MICRO) {
+-                            tiled_run_microtile<SUBBLK, HAS_MIN, BIAS, ACTBIAS>(ws->src0, ws->src1,
++                            tiled_run_microtile<SUBBLK, HAS_MIN, BIAS, ACTBIAS, Q41>(ws->src0, ws->src1,
+                                 (int) (ir0 - iir0), (int) (ir1 - iir1), 1, 0,
+                                 ws->acc, TILED_TILE_ROWS);
+                         }
+@@ -1024,7 +1086,7 @@ static void ggml_compute_forward_mul_mat_tiled_one_chunk(
+     }
+ }
+ 
+-template <typename B, int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS>
++template <typename B, int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS, typename BA = block_q8_K, bool Q41 = false>
+ static void ggml_compute_forward_mul_mat_tiled_driver(
+         const struct ggml_compute_params * params,
+               struct ggml_tensor * dst) {
+@@ -1133,7 +1195,7 @@ static void ggml_compute_forward_mul_mat_tiled_driver(
+         const int64_t ir1_start = dr1 * ith1;
+         const int64_t ir1_end = MIN(ir1_start + dr1, nr1);
+ 
+-        ggml_compute_forward_mul_mat_tiled_one_chunk<B, SUBBLK, HAS_MIN, BIAS, ACTBIAS>(params, dst, ir0_start, ir0_end, ir1_start, ir1_end, ws);
++        ggml_compute_forward_mul_mat_tiled_one_chunk<B, SUBBLK, HAS_MIN, BIAS, ACTBIAS, BA, Q41>(params, dst, ir0_start, ir0_end, ir1_start, ir1_end, ws);
+ 
+         if (nth >= nchunk0 * nchunk1) {
+             break;
+@@ -1145,7 +1207,7 @@ static void ggml_compute_forward_mul_mat_tiled_driver(
+ 
+ // src0 type dispatch, shared by the MUL_MAT and MUL_MAT_ID entries: one expert for
+ // MUL_MAT_ID (expert_rows != NULL), the full op for MUL_MAT
+-template <typename B, int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS>
++template <typename B, int SUBBLK, bool HAS_MIN, int BIAS, bool ACTBIAS, typename BA = block_q8_K, bool Q41 = false>
+ static bool tiled_matmul_dispatch(const struct ggml_compute_params * params,
+                                   struct ggml_tensor * dst,
+                                   const int32_t * expert_rows,
+@@ -1153,9 +1215,9 @@ static bool tiled_matmul_dispatch(const struct ggml_compute_params * params,
+                                   int64_t cne1,
+                                   char * scratch) {
+     if (expert_rows == NULL) {
+-        ggml_compute_forward_mul_mat_tiled_driver<B, SUBBLK, HAS_MIN, BIAS, ACTBIAS>(params, dst);
++        ggml_compute_forward_mul_mat_tiled_driver<B, SUBBLK, HAS_MIN, BIAS, ACTBIAS, BA, Q41>(params, dst);
+     } else {
+-        ggml_compute_forward_mul_mat_id_tiled_one_expert<B, SUBBLK, HAS_MIN, BIAS, ACTBIAS>(params, dst, cur_a, cne1, expert_rows, scratch);
++        ggml_compute_forward_mul_mat_id_tiled_one_expert<B, SUBBLK, HAS_MIN, BIAS, ACTBIAS, BA, Q41>(params, dst, cur_a, cne1, expert_rows, scratch);
+     }
+     return true;
+ }
+@@ -1174,6 +1236,8 @@ static bool ggml_tiled_matmul_type_dispatch(const struct ggml_compute_params * p
+             return tiled_matmul_dispatch<block_q5_K, 32, true,  0, false>(params, dst, expert_rows, cur_a, cne1, scratch);
+         case GGML_TYPE_Q4_K:
+             return tiled_matmul_dispatch<block_q4_K, 32, true,  0, false>(params, dst, expert_rows, cur_a, cne1, scratch);
++        case GGML_TYPE_Q4_1:
++            return tiled_matmul_dispatch<block_q4_1, 32, false, 0, false, block_q8_1, true>(params, dst, expert_rows, cur_a, cne1, scratch);
+         case GGML_TYPE_Q3_K:
+             return tiled_matmul_dispatch<block_q3_K, 16, false,  4, true>(params, dst, expert_rows, cur_a, cne1, scratch);
+         case GGML_TYPE_Q2_K:
 diff --git a/ggml/src/ggml.c b/ggml/src/ggml.c
 index 5d8449a88..203f15e43 100644
 --- a/ggml/src/ggml.c