q4_K mat-vec tuning for gfx1151

ggml/src/ggml-cuda/common.cuh

@@ -371,6 +371,16 @@ static bool ggml_cuda_is_aligned(const ggml_tensor * tensor, const size_t alignm
            tensor->nb[3] % alignment == 0;
 }
 
+// q8_1 activations packed 4 blocks at a time, mirroring ggml-vulkan's block_q8_1_x4.
+// The scales are hoisted out of the per-block headers so that the quant bytes are
+// contiguous and 16 B aligned, which allows a 128-bit activation load. Same total size
+// and same byte offset for any 4-block-aligned index as an array of block_q8_1.
+struct block_q8_1_x4 {
+    half2  ds[4];    //  16 B
+    int8_t qs[128];  // 128 B: four blocks' quant arrays, concatenated
+};
+static_assert(sizeof(block_q8_1_x4) == 4*sizeof(block_q8_1), "block_q8_1_x4 must alias 4 q8_1 blocks");
+
 static constexpr __device__ int ggml_cuda_get_physical_warp_size() {
 #if defined(GGML_USE_HIP) && (defined(__GFX9__) || defined(__GFX8__))
     return 64;

ggml/src/ggml-cuda/mmvq.cu

@@ -349,7 +349,25 @@ static constexpr __device__ int get_mmvq_mmid_max_batch_for_device() {
 #endif
 }
 
+// Manually unroll the q4_K K-loop, mirroring ggml-vulkan's mul_mat_vecq.comp. The trip count
+// is runtime-dependent so the compiler will not unroll it on its own; doing it by hand puts
+// four iterations of loads+dots in straight-line code, raising memory-level parallelism.
+// The sched_group_barrier is not optional: unrolling without it is markedly slower than not
+// unrolling at all, because the whole unrolled body's live set is scheduled at once.
+#if defined(GGML_USE_HIP) && defined(RDNA3_5)
+#define GGML_MMVQ_Q4K_UNROLL 4
+#define GGML_MMVQ_Q4K_SCHED_GROUP_BARRIER()                    \
+    do {                                                       \
+        __builtin_amdgcn_sched_group_barrier(0x020, 2, 0);     \
+        __builtin_amdgcn_sched_group_barrier(0x002, 8, 0);     \
+    } while (0)
+#endif
+
 static constexpr __host__ __device__ int calc_nwarps(ggml_type type, int ncols_dst, mmvq_parameter_table_id table_id) {
+    if (table_id == MMVQ_PARAMETERS_RDNA2 && ncols_dst == 1 && type == GGML_TYPE_Q4_K) {
+        // one wave per output row leaves only 2 q4_K superblocks in flight; 2 measures faster
+        return 2;
+    }
     if (table_id == MMVQ_PARAMETERS_GENERIC) {
         switch (ncols_dst) {
             case 1:
@@ -589,11 +607,11 @@ static __global__ void mul_mat_vec_q(
     const block_q8_1 * y = ((const block_q8_1 *) vy) + sample_y*stride_sample_y + channel_y*stride_channel_y;
     const int kbx_offset = sample_x*stride_sample_x + channel_x*stride_channel_x + row0*stride_row_x;
 
-    for (int kbx = tid / (qi/vdr); kbx < blocks_per_row_x; kbx += blocks_per_iter) {
-        const int kby = kbx * (qk/QK8_1); // y block index that aligns with kbx
+    // x block quant index when casting the quants to int
+    const int kqs = vdr * (tid % (qi/vdr));
 
-        // x block quant index when casting the quants to int
-        const int kqs = vdr * (tid % (qi/vdr));
+    auto iter = [&](const int kbx) {
+        const int kby = kbx * (qk/QK8_1); // y block index that aligns with kbx
 
 #pragma unroll
         for (int j = 0; j < ncols_dst; ++j) {
@@ -609,6 +627,33 @@ static __global__ void mul_mat_vec_q(
                 }
             }
         }
+    };
+
+#ifdef GGML_MMVQ_Q4K_UNROLL
+    if constexpr (type == GGML_TYPE_Q4_K) {
+        const int kbx0 = tid / (qi/vdr);
+        int kbx  = kbx0;
+        int n_it = kbx0 < blocks_per_row_x
+                 ? (blocks_per_row_x - kbx0 + blocks_per_iter - 1) / blocks_per_iter : 0;
+        while (n_it >= GGML_MMVQ_Q4K_UNROLL) {
+#pragma unroll
+            for (int u = 0; u < GGML_MMVQ_Q4K_UNROLL; ++u) {
+                iter(kbx);
+                kbx += blocks_per_iter;
+                GGML_MMVQ_Q4K_SCHED_GROUP_BARRIER();
+            }
+            n_it -= GGML_MMVQ_Q4K_UNROLL;
+        }
+        while (n_it-- > 0) {
+            iter(kbx);
+            kbx += blocks_per_iter;
+        }
+    } else
+#endif // GGML_MMVQ_Q4K_UNROLL
+    {
+        for (int kbx = tid / (qi/vdr); kbx < blocks_per_row_x; kbx += blocks_per_iter) {
+            iter(kbx);
+        }
     }
 
     __shared__ float tmp_shared[nwarps-1 > 0 ? nwarps-1 : 1][ncols_dst][rows_per_cuda_block][warp_size];

ggml/src/ggml-cuda/quantize.cu

@@ -50,6 +50,7 @@ static __device__ __forceinline__ float nvfp4_native_scale_error(
 }
 #endif // CUDART_VERSION >= 12080
 
+template <bool use_x4>
 __launch_bounds__(CUDA_QUANTIZE_BLOCK_SIZE, 1)
 static __global__ void quantize_q8_1(
         const float * x_ptr, void * vy_ptr,
@@ -91,13 +92,36 @@ static __global__ void quantize_q8_1(
     const float  d = amax / 127.0f;
     const int8_t q = amax == 0.0f ? 0 : roundf(xi / d);
 
-    y[ib].qs[iqs] = q;
+    if constexpr (use_x4) {
+        // Scales hoisted to the front of each 4-block group, quant bytes contiguous.
+        //
+        // ds.y holds d*sum(q), the DEQUANTIZED block sum, which is the convention
+        // ggml-vulkan's quantize_q8_1 uses and what the q4_K min term below consumes. The
+        // plain block_q8_1 layout stores sum(x) instead; the two differ by the quantization
+        // residual, and feeding sum(x) to that min term costs an order of magnitude of
+        // accuracy on long rows.
+        const float sumq = warp_reduce_sum<QK8_1>((float) q);
 
-    if (iqs > 0) {
-        return;
-    }
+        block_q8_1_x4 * y4 = (block_q8_1_x4 *) vy;
+        const int64_t outer = ib >> 2;
+        const int64_t inner = ib &  3;
+
+        y4[outer].qs[inner*QK8_1 + iqs] = q;
+
+        if (iqs > 0) {
+            return;
+        }
+
+        y4[outer].ds[inner] = make_half2(d, d*sumq);
+    } else {
+        y[ib].qs[iqs] = q;
+
+        if (iqs > 0) {
+            return;
+        }
 
-    y[ib].ds = make_half2(d, sum);
+        y[ib].ds = make_half2(d, sum);
+    }
 }
 
 __device__ __forceinline__ uint8_t compute_e8m0_scale(float amax) {
@@ -568,8 +592,14 @@ void quantize_row_q8_1_cuda(
     const dim3 num_blocks(block_num_x, ne1, ne2*ne3);
     const dim3 block_size(CUDA_QUANTIZE_BLOCK_SIZE, 1, 1);
     const ggml_cuda_kernel_launch_params launch_params = ggml_cuda_kernel_launch_params(num_blocks, block_size, 0, stream);
-    ggml_cuda_kernel_launch(quantize_q8_1, launch_params, x, vy, ne00, s01, s02, s03, ne0, ne1, ne2_fastdiv);
-    GGML_UNUSED(type_src0);
+    // Only the RDNA3.5 q4_K vec_dot reads the x4 layout; every other type and arch keeps the
+    // plain block_q8_1 layout. The buffer is allocated per mul_mat, so the two can coexist.
+    const int cc = ggml_cuda_info().devices[ggml_cuda_get_device()].cc;
+    if (type_src0 == GGML_TYPE_Q4_K && GGML_CUDA_CC_IS_RDNA3_5(cc)) {
+        ggml_cuda_kernel_launch(quantize_q8_1<true>,  launch_params, x, vy, ne00, s01, s02, s03, ne0, ne1, ne2_fastdiv);
+    } else {
+        ggml_cuda_kernel_launch(quantize_q8_1<false>, launch_params, x, vy, ne00, s01, s02, s03, ne0, ne1, ne2_fastdiv);
+    }
 }
 
 void quantize_mmq_q8_1_cuda(

ggml/src/ggml-cuda/vecdotq.cuh

@@ -891,6 +891,72 @@ static __device__ __forceinline__ float vec_dot_q4_K_q8_1(
 
     const block_q4_K * bq4_K = (const block_q4_K *) vbq + kbx;
 
+#if defined(RDNA3_5)
+    // Mirrors ggml-vulkan's mul_mat_vecq.comp: each thread takes one aligned 16-byte chunk
+    // of qs and one nibble half, giving a 128-bit weight load, plus 16 contiguous activation
+    // bytes for a 128-bit activation load (requires the block_q8_1_x4 layout). The two threads
+    // sharing a chunk read the same 16 bytes with different shifts, the same register-level
+    // redundancy RADV accepts in exchange for wide loads. Still 16 threads/superblock, so VDR
+    // stays 2 and the K-loop trip count is unchanged.
+    const int j = iqs >> 1;            // 0..15
+    const int c = j >> 1;              // 16-byte chunk of qs
+    const int h = j &  1;              // nibble half
+
+    // 16 contiguous weights starting at W, all inside one 32-weight sub-block
+    const int W  = 64*(c >> 1) + 16*(c & 1) + 32*h;
+    const int sb = W >> 5;             // sub-block / q8_1 block index, 0..7
+    const int wo = W & 31;             // byte offset inside that block, 0 or 16
+
+    // qs is at offset 16 in a 144-byte block, so qs + 16*c is always 16 B aligned
+    const uint4 wv = *(const uint4 *) __builtin_assume_aligned(bq4_K->qs + 16*c, 16);
+    const int sh = 4*h;
+    const int w0 = (wv.x >> sh) & 0x0F0F0F0F;
+    const int w1 = (wv.y >> sh) & 0x0F0F0F0F;
+    const int w2 = (wv.z >> sh) & 0x0F0F0F0F;
+    const int w3 = (wv.w >> sh) & 0x0F0F0F0F;
+
+    const block_q8_1_x4 * bq8x4 = (const block_q8_1_x4 *) bq8_1;
+    const int8_t * qs8 = bq8x4[sb >> 2].qs + (sb & 3)*QK8_1 + wo;
+    const uint4 uv = *(const uint4 *) __builtin_assume_aligned(qs8, 16);
+
+    const int u0 = uv.x, u1 = uv.y, u2 = uv.z, u3 = uv.w;
+
+    int sumi_d = 0;
+    sumi_d = ggml_cuda_dp4a(w0, u0, sumi_d);
+    sumi_d = ggml_cuda_dp4a(w1, u1, sumi_d);
+    sumi_d = ggml_cuda_dp4a(w2, u2, sumi_d);
+    sumi_d = ggml_cuda_dp4a(w3, u3, sumi_d);
+
+    // Branchless get_scale_min_k4 for sub-block sb. Only three distinct bytes are ever needed;
+    // load them unconditionally and select, which keeps this in VALU instead of emitting
+    // branches (or conditional loads) in the hot loop.
+    const uint8_t * sc8 = bq4_K->scales;
+    const int hi = sb >> 2;
+    const int A = sc8[sb];
+    const int B = sc8[sb + 4];
+    const int C = sc8[sb & 3];
+    const int s_a = hi ? B : A;
+    const int s_b = B;
+    const int s_c = hi ? C : A;
+    const int s_s = A;
+    const int sc_lo = s_s & 63;
+    const int mn_lo = s_b & 63;
+    const int sc_hi = (s_a & 0x0F) | ((s_c >> 6) << 4);
+    const int mn_hi = (s_a >>   4) | ((s_s >> 6) << 4);
+    // NOTE: keep these signed. sumi_d is a dot product of signed int8 and is frequently
+    // negative; an unsigned scale would make sumi_d*sc unsigned arithmetic.
+    const int sc = hi ? sc_hi : sc_lo;
+    const int mn = hi ? mn_hi : mn_lo;
+
+    // The min term needs sum(u) over this thread's 16 activations. q8_1 already carries the
+    // whole 32-element block sum in ds.y, so read it instead of recomputing it with 4 more
+    // dp4a. Per lane the ds.y*0.5 split is approximate, but the two threads sharing a q8_1
+    // block also share the sub-block scale, so it is exact after the cross-lane reduction.
+    const float2 ds8  = __half22float2(bq8x4[sb >> 2].ds[sb & 3]);
+    const float2 dm4f = __half22float2(bq4_K->dm);
+    return dm4f.x * (ds8.x * (sumi_d * sc)) - dm4f.y * (mn * ds8.y * 0.5f);
+#else
+
     int    v[2];
     int    u[2*QR4_K];
     float d8[QR4_K];
@@ -930,6 +996,7 @@ static __device__ __forceinline__ float vec_dot_q4_K_q8_1(
     }
 
     return vec_dot_q4_K_q8_1_impl_vmmq(v, u, sc, m, bq4_K->dm, d8);
+#endif // defined(RDNA3_5)
 }
 
 static __device__ __forceinline__ float vec_dot_q5_K_q8_1(

tests/test-backend-ops.cpp

@@ -8069,6 +8069,13 @@ static void add_rdna35_mmq_cases(std::vector<std::unique_ptr<test_case>> & test_
 
     test_cases.emplace_back(new test_mul_mat(GGML_TYPE_Q4_K, GGML_TYPE_F32,  512, 16, 2048, {1, 1}, {1, 1}));
     test_cases.emplace_back(new test_mul_mat(GGML_TYPE_Q4_K, GGML_TYPE_F32, 4096, 16, 4096, {1, 1}, {1, 1}));
+
+    // gemma-4-31B-it Q4_K_M decode: the FFN gate/up matvec, the largest single op in decode.
+    // In the model this op is fused (MUL_MAT + MUL_MAT + GLU), which a plain test_mul_mat does
+    // not reproduce; this case exists to time the mul_mat_vec_q kernel itself, for which the
+    // perf harness duplicates the MUL_MAT node directly.
+    test_cases.emplace_back(new test_mul_mat(GGML_TYPE_Q4_K, GGML_TYPE_F32, 21504, 1, 5376, {1, 1}, {1, 1}));
+
     test_cases.emplace_back(new test_mul_mat(GGML_TYPE_Q6_K, GGML_TYPE_F32, 4096, 128, 12288, {1, 1}, {1, 1}));
     test_cases.emplace_back(new test_mul_mat(GGML_TYPE_Q8_0, GGML_TYPE_F32, 32, 128, 4096, {1, 1}, {1, 1}));
     test_cases.emplace_back(new test_mul_mat(GGML_TYPE_Q4_0, GGML_TYPE_F32, 32, 128, 4096, {1, 1}, {1, 1}));