tests : judge the batched DT3 path strictly when it is MMQ

With MMQ, ncols_dst > 8 is an integer path in the same numerical regime
as MMVQ and is judged against the exact reference at 1e-5 instead of
riding the loose GEMM gate. The regime is told apart by the result
itself; a GEMM fallback (backends without DT3 MMQ) keeps the fp16
reference and the F16 GEMM bit-identity control. Also add ncols_dst=100
to exercise a wide tile with a clamped last column block.

ggml/src/ggml-cuda/mmq-config-ampere.cuh

@@ -33,6 +33,23 @@ static constexpr __host__ __device__ ggml_cuda_mmq_config ggml_cuda_mmq_get_conf
     CASE(GGML_TYPE_Q2_0, 256, 1, 128, 112, GGML_CUDA_MMQ_SRAM_LAYOUT_Q8_0, MMQ_ITER_K, true, false);
     CASE(GGML_TYPE_Q2_0, 256, 1, 128, 128, GGML_CUDA_MMQ_SRAM_LAYOUT_Q8_0, MMQ_ITER_K, true, false);
 
+    CASE(GGML_TYPE_DT3, 256, 1, 128,   8, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, true);
+    CASE(GGML_TYPE_DT3, 256, 1, 128,  16, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, true);
+    CASE(GGML_TYPE_DT3, 256, 1, 128,  32, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, true);
+    CASE(GGML_TYPE_DT3, 256, 1, 128,  64, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, true);
+    CASE(GGML_TYPE_DT3, 256, 1, 128, 128, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, true);
+    CASE(GGML_TYPE_DT3, 256, 1, 128,   8, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, false);
+    CASE(GGML_TYPE_DT3, 256, 1, 128,  16, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, false);
+    CASE(GGML_TYPE_DT3, 256, 1, 128,  24, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, false);
+    CASE(GGML_TYPE_DT3, 256, 1, 128,  32, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, false);
+    CASE(GGML_TYPE_DT3, 256, 1, 128,  40, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, false);
+    CASE(GGML_TYPE_DT3, 256, 1, 128,  48, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, false);
+    CASE(GGML_TYPE_DT3, 256, 1, 128,  64, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, false);
+    CASE(GGML_TYPE_DT3, 256, 1, 128,  80, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, false);
+    CASE(GGML_TYPE_DT3, 256, 1, 128,  96, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, false);
+    CASE(GGML_TYPE_DT3, 256, 1, 128, 112, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, false);
+    CASE(GGML_TYPE_DT3, 256, 1, 128, 128, GGML_CUDA_MMQ_SRAM_LAYOUT_DT3, MMQ_ITER_K, true, false);
+
     CASE(GGML_TYPE_Q4_0, 256, 1, 128,   8, GGML_CUDA_MMQ_SRAM_LAYOUT_Q8_0, MMQ_ITER_K, true, true);
     CASE(GGML_TYPE_Q4_0, 256, 1, 128,  16, GGML_CUDA_MMQ_SRAM_LAYOUT_Q8_0, MMQ_ITER_K, true, true);
     CASE(GGML_TYPE_Q4_0, 256, 1, 128,  32, GGML_CUDA_MMQ_SRAM_LAYOUT_Q8_0, MMQ_ITER_K, true, true);

ggml/src/ggml-cuda/mmq-load-tiles.cuh

@@ -176,6 +176,117 @@ template <ggml_type type, int J, bool fallback> static __device__ __forceinline_
     }
 }
 
+template <ggml_type type, int J, bool fallback> static __device__ __forceinline__ void ggml_cuda_mmq_load_tiles_dt3(
+        const char * __restrict__ x, int * __restrict__ x_tile, const int kbx0, const int i_max, const int stride) {
+    constexpr int warp_size   = ggml_cuda_get_physical_warp_size();
+    constexpr int nwarps      = ggml_cuda_mmq_get_nthreads(type, J, fallback) / warp_size;
+    constexpr int I           = ggml_cuda_mmq_get_I(type, J, fallback);
+    constexpr int sram_stride = ggml_cuda_mmq_get_sram_stride(type, J, fallback);
+
+    // DT3: 128 elements as two ternary planes with one fp16 scale each, w = d1*t1 + d2*t2.
+    // The two planes cannot be fused into a single int8 value because d1 != d2, so the
+    // tile holds both planes decoded to int8 trits in {-1, 0, +1}, plane 2 offset by
+    // 2*MMQ_TILE_NE_K ints from plane 1 within each row, and 2x8 float scales per row.
+    // The decode is the same base-3 digit iteration as vec_dot_dt3_q8_1 (see vecdotq.cuh):
+    // each packed byte is decoded once with q -> (q*3) & 0xFF, two bytes at a time in the
+    // 16-bit lanes of one int. Digit bytes in {0, 1, 2} become trit bytes in {-1, 0, +1}
+    // without cross-byte borrows via ((dig | 0x80808080) - 0x01010101) ^ 0x80808080.
+#if defined(AMD_MFMA_AVAILABLE) || defined(TURING_MMA_AVAILABLE) || defined(AMD_WMMA_AVAILABLE)
+    constexpr int row_stride_qs = sram_stride;
+    constexpr int row_stride_df = sram_stride;
+    int   * x_qs = (int   *)  x_tile;
+    float * x_df = (float *) (x_qs + 4*MMQ_TILE_NE_K);
+#else
+    constexpr tile_x_sizes txs = mmq_get_dp4a_tile_x_sizes(GGML_TYPE_DT3, I);
+    constexpr int row_stride_qs = 4*MMQ_TILE_NE_K + 1;
+    constexpr int row_stride_df = 4*MMQ_TILE_NE_K/QI8_0 + 1;
+    int   * x_qs = (int   *)  x_tile;
+    float * x_df = (float *) (x_qs + txs.qs);
+#endif // defined(AMD_MFMA_AVAILABLE) || defined(TURING_MMA_AVAILABLE) || defined(AMD_WMMA_AVAILABLE)
+
+    constexpr int blocks_per_iter = MMQ_ITER_K / QK_DT3;
+    static_assert(blocks_per_iter == 2, "DT3 load assumes 2 blocks per iteration");
+
+    // 32 threads per row: 2 blocks x 2 planes x 8 slots. Slots 0..5 decode one 4-byte
+    // quad of qs each (5 ints of 4 trits), slot 6 decodes the 2 qh bytes (2 ints),
+    // slot 7 is idle.
+    constexpr int threads_per_row = 32;
+    constexpr int nrows = warp_size / threads_per_row;
+
+    const int txi = threadIdx.x % threads_per_row;
+    const int kbx = txi / 16;
+    const int p   = (txi / 8) % 2;
+    const int u   = txi % 8;
+
+#pragma unroll
+    for (int i0 = 0; i0 < I; i0 += nrows*nwarps) {
+        int i = i0 + threadIdx.y*nrows + threadIdx.x/threads_per_row;
+
+        if (fallback) {
+            i = min(i, i_max);
+        }
+
+        const block_dt3 * bxi = (const block_dt3 *) x + kbx0 + i*stride + kbx;
+        int * dst = x_qs + i*row_stride_qs + p*(2*MMQ_TILE_NE_K) + kbx*(QK_DT3/4);
+
+        if (u < 6) {
+            // qs[0..16): byte m holds elements m + 16*n; qs[16..24): byte 16 + m
+            // holds elements 80 + m + 8*n. Either way a quad of consecutive bytes
+            // yields 4 consecutive elements per digit n, i.e. one tile int.
+            const int q32 = get_int_b4(bxi->qs[p], u);
+            int qa = (q32 >> 0) & 0x00FF00FF;
+            int qb = (q32 >> 8) & 0x00FF00FF;
+#pragma unroll
+            for (int n = 0; n < 5; ++n) {
+                const int qa3 = qa*3;
+                const int qb3 = qb*3;
+                const int dig = ((qa3 >> 8) & 0x00030003) | (qb3 & 0x03000300);
+                qa = qa3 & 0x00FF00FF;
+                qb = qb3 & 0x00FF00FF;
+
+                const int idx = u < 4 ? 4*n + u : 20 + 2*n + (u - 4);
+                dst[idx] = ((dig | 0x80808080) - 0x01010101) ^ 0x80808080;
+            }
+        } else if (u == 6) {
+            // qh: byte b holds elements 120 + b + 2*n for n = 0..3 — only 4 digits
+            // are iterated, the 5th is packing padding and must never decode.
+            int q = bxi->qh[p][0] | (bxi->qh[p][1] << 16);
+#pragma unroll
+            for (int s = 0; s < 2; ++s) {
+                int dig = 0;
+#pragma unroll
+                for (int n = 0; n < 2; ++n) {
+                    const int q3 = q*3;
+                    dig |= (((q3 >> 8) & 0x03) | ((q3 >> 16) & 0x0300)) << (16*n);
+                    q = q3 & 0x00FF00FF;
+                }
+                dst[30 + s] = ((dig | 0x80808080) - 0x01010101) ^ 0x80808080;
+            }
+        }
+    }
+
+    // 16 scale entries per row and iteration: 2 planes x 2 blocks x 4 q8_1 chunks.
+    // Plane 2 scales sit after the 8 plane-1 entries, matching the vec_dot indexing.
+    constexpr int scale_entries_per_plane = blocks_per_iter*(QK_DT3/QK8_1);
+
+    const int ksx = threadIdx.x % (2*scale_entries_per_plane);
+    const int ps  = ksx / scale_entries_per_plane;
+    const int scale_block = (ksx % scale_entries_per_plane) / (QK_DT3/QK8_1);
+
+#pragma unroll
+    for (int i0 = 0; i0 < I; i0 += nwarps) {
+        int i = i0 + threadIdx.y;
+
+        if (fallback) {
+            i = min(i, i_max);
+        }
+
+        const block_dt3 * bxi = (const block_dt3 *) x + kbx0 + i*stride + scale_block;
+
+        x_df[i*row_stride_df + ksx] = bxi->d[ps];
+    }
+}
+
 template <ggml_type type, int J, bool fallback> static __device__ __forceinline__ void ggml_cuda_mmq_load_tiles_q4_0(
         const char * __restrict__ x, int * __restrict__ x_tile, const int kbx0, const int i_max, const int stride) {
     constexpr int warp_size   = ggml_cuda_get_physical_warp_size();

ggml/src/ggml-cuda/mmq-vec-dot.cuh

@@ -280,6 +280,140 @@ static __device__ __forceinline__ void ggml_cuda_mmq_vec_dot_q8_0_q8_1_mma(
 }
 
 
+// DT3: both decoded ternary planes of the tile are multiplied against the same y data,
+// each with its own per-chunk scale: sum = dB * (sumi1*dA1 + sumi2*dA2).
+template <ggml_type type, int J, bool fallback> static __device__ __forceinline__ void ggml_cuda_mmq_vec_dot_dt3_q8_1_dp4a(
+        const int * __restrict__ x, const int * __restrict__ y, float * __restrict__ sum, const int k00) {
+    constexpr int warp_size = ggml_cuda_get_physical_warp_size();
+    constexpr int nwarps    = ggml_cuda_mmq_get_nthreads(type, J, fallback) / warp_size;
+    constexpr int I         = ggml_cuda_mmq_get_I(type, J, fallback);
+
+    constexpr tile_x_sizes txs = mmq_get_dp4a_tile_x_sizes(GGML_TYPE_DT3, I);
+    const int   * x_qs = (const int   *) x;
+    const float * x_df = (const float *) x_qs + txs.qs;
+    const int   * y_qs = (const int   *) y + 4;
+    const float * y_df = (const float *) y;
+
+    constexpr int row_stride_qs = 4*MMQ_TILE_NE_K + 1;
+    constexpr int row_stride_df = 4*MMQ_TILE_NE_K/QI8_0 + 1;
+
+// #pragma unroll
+    for (int k01 = 0; k01 < MMQ_TILE_NE_K; k01 += VDR_Q8_0_Q8_1_MMQ) {
+        const int k0 = k00 + k01;
+
+#pragma unroll
+        for (int j0 = 0; j0 < J; j0 += nwarps) {
+            const int j = j0 + threadIdx.y;
+
+#pragma unroll
+            for (int i0 = 0; i0 < I; i0 += warp_size) {
+                const int i = i0 + threadIdx.x;
+
+                const int   * yqs = &y_qs[j*MMQ_TILE_Y_K + k0 % MMQ_TILE_NE_K];
+                const float   dB  =  y_df[j*MMQ_TILE_Y_K + (k0/QI8_1) % (MMQ_TILE_NE_K/QI8_1)];
+
+                sum[j0/nwarps*I/warp_size + i0/warp_size] += vec_dot_q8_0_q8_1_impl<float, VDR_Q8_0_Q8_1_MMQ>
+                    (&x_qs[i*row_stride_qs + k0], yqs,
+                     x_df[i*row_stride_df + k0/QI8_0], dB);
+                sum[j0/nwarps*I/warp_size + i0/warp_size] += vec_dot_q8_0_q8_1_impl<float, VDR_Q8_0_Q8_1_MMQ>
+                    (&x_qs[i*row_stride_qs + 2*MMQ_TILE_NE_K + k0], yqs,
+                     x_df[i*row_stride_df + 2*MMQ_TILE_NE_K/QI8_0 + k0/QI8_0], dB);
+            }
+        }
+    }
+}
+
+template <ggml_type type, int J, bool fallback>
+static __device__ __forceinline__ void ggml_cuda_mmq_vec_dot_dt3_q8_1_mma(
+        const int * __restrict__ x, const int * __restrict__ y, float * __restrict__ sum, const int k00) {
+#if defined(AMD_MFMA_AVAILABLE) || defined(AMD_WMMA_AVAILABLE)
+    // DT3 MMQ is not enabled for AMD (no config entries select it); this only has to compile.
+    GGML_UNUSED_VARS(x, y, sum, k00);
+    NO_DEVICE_CODE;
+#else
+    typedef tile<16, 8, int> tile_A;
+    typedef tile< 8, 8, int> tile_B;
+    typedef tile<16, 8, int> tile_C;
+
+    constexpr int I             = ggml_cuda_mmq_get_I(type, J, fallback);
+    constexpr int sram_stride   = ggml_cuda_mmq_get_sram_stride(type, J, fallback);
+    constexpr int rows_per_warp = ggml_cuda_mmq_get_rows_per_warp(type, J, fallback);
+    constexpr int ntx           = rows_per_warp/tile_C::I; // Number of x minitiles per warp.
+
+    y += (threadIdx.y % ntx) * (tile_C::J*MMQ_TILE_Y_K);
+
+    const int   * x_qs = (const int   *) x;
+    const float * x_df = (const float *) x_qs + 4*MMQ_TILE_NE_K;
+    const int   * y_qs = (const int   *) y + 4;
+    const float * y_df = (const float *) y;
+
+    tile_A A[ntx][2][MMQ_TILE_NE_K/QI8_0];
+    float dA[ntx][tile_C::ne/2][2][MMQ_TILE_NE_K/QI8_0];
+
+    const int i0 = (threadIdx.y/ntx)*rows_per_warp;
+
+#pragma unroll
+    for (int n = 0; n < ntx; ++n) {
+#pragma unroll
+        for (int p = 0; p < 2; ++p) {
+#pragma unroll
+            for (int k01 = 0; k01 < MMQ_TILE_NE_K; k01 += QI8_0) {
+                const int k0 = k00 + k01;
+
+                load_ldmatrix(A[n][p][k01/QI8_0], x_qs + (i0 + n*tile_A::I)*sram_stride + p*(2*MMQ_TILE_NE_K) + k0, sram_stride);
+            }
+        }
+
+#pragma unroll
+        for (int l = 0; l < tile_C::ne/2; ++l) {
+            const int i = i0 + n*tile_A::I + tile_C::get_i(2*l);
+
+#pragma unroll
+            for (int p = 0; p < 2; ++p) {
+#pragma unroll
+                for (int k01 = 0; k01 < MMQ_TILE_NE_K; k01 += QI8_0) {
+                    const int k0 = k00 + k01;
+
+                    dA[n][l][p][k01/QI8_0] = x_df[i*sram_stride + p*(2*MMQ_TILE_NE_K/QI8_0) + k0/QI8_0];
+                }
+            }
+        }
+    }
+
+#pragma unroll
+    for (int j0 = 0; j0 < J; j0 += ntx*tile_C::J) {
+#pragma unroll
+        for (int k01 = 0; k01 < MMQ_TILE_NE_K; k01 += QI8_0) {
+            tile_B B;
+            float dB[tile_C::ne/2];
+
+            load_generic(B, y_qs + j0*MMQ_TILE_Y_K + k01, MMQ_TILE_Y_K); // faster than load_ldmatrix
+
+#pragma unroll
+            for (int l = 0; l < tile_C::ne/2; ++l) {
+                const int j = j0 + tile_C::get_j(l);
+
+                dB[l] = y_df[j*MMQ_TILE_Y_K + k01/QI8_1];
+            }
+
+#pragma unroll
+            for (int n = 0; n < ntx; ++n) {
+                tile_C C1;
+                tile_C C2;
+                mma(C1, A[n][0][k01/QI8_0], B);
+                mma(C2, A[n][1][k01/QI8_0], B);
+
+#pragma unroll
+                for (int l = 0; l < tile_C::ne; ++l) {
+                    sum[(j0/tile_C::J + n)*tile_C::ne + l] +=
+                        (C1.x[l]*dA[n][l/2][0][k01/QI8_0] + C2.x[l]*dA[n][l/2][1][k01/QI8_0])*dB[l%2];
+                }
+            }
+        }
+    }
+#endif // defined(AMD_MFMA_AVAILABLE) || defined(AMD_WMMA_AVAILABLE)
+}
+
 template <ggml_type type, int J, bool fallback> static __device__ __forceinline__ void ggml_cuda_mmq_vec_dot_q8_1_q8_1_dp4a(
         const int * __restrict__ x, const int * __restrict__ y, float * __restrict__ sum, const int k00) {
     constexpr int warp_size = ggml_cuda_get_physical_warp_size();

ggml/src/ggml-cuda/mmq.cu

@@ -13,6 +13,9 @@ static void ggml_cuda_mul_mat_q_switch_type(ggml_backend_cuda_context & ctx, con
         case GGML_TYPE_Q2_0:
             mul_mat_q_case<GGML_TYPE_Q2_0>(ctx, args, stream);
             break;
+        case GGML_TYPE_DT3:
+            mul_mat_q_case<GGML_TYPE_DT3>(ctx, args, stream);
+            break;
         case GGML_TYPE_Q4_0:
             mul_mat_q_case<GGML_TYPE_Q4_0>(ctx, args, stream);
             break;
@@ -261,6 +264,18 @@ bool ggml_cuda_should_use_mmq(enum ggml_type type, int cc, int64_t ne11, int64_t
     return false;
 #endif // GGML_CUDA_FORCE_CUBLAS
 
+    // DT3 keeps two decoded ternary planes per row in SRAM, roughly double the tile
+    // of a single-plane type: only the MMA data layout is implemented and even the
+    // narrowest tile needs ~76 KiB of shared memory, more than e.g. Turing offers.
+    if (type == GGML_TYPE_DT3) {
+        if (!turing_mma_available(cc)) {
+            return false;
+        }
+        const int    id    = ggml_cuda_get_device();
+        const size_t smpbo = ggml_cuda_info().devices[id].smpbo;
+        return mmq_get_nbytes_shared(ggml_cuda_mmq_get_config(GGML_TYPE_DT3, 8, true, cc), cc) <= smpbo;
+    }
+
     bool mmq_supported;
 
     switch (type) {

ggml/src/ggml-cuda/mmq.cuh

@@ -61,6 +61,7 @@ static mmq_q8_1_ds_layout mmq_get_q8_1_ds_layout(const ggml_type type_x) {
     switch (type_x) {
         case GGML_TYPE_Q1_0:
         case GGML_TYPE_Q2_0:
+        case GGML_TYPE_DT3:
             return MMQ_Q8_1_DS_LAYOUT_D4;
         case GGML_TYPE_Q4_0:
         case GGML_TYPE_Q4_1:
@@ -121,6 +122,7 @@ struct tile_x_sizes {
 
 enum ggml_cuda_mmq_sram_layout {
     GGML_CUDA_MMQ_SRAM_LAYOUT_Q8_0,
+    GGML_CUDA_MMQ_SRAM_LAYOUT_DT3,   // Two decoded ternary planes per row, each with its own per-block scales.
     GGML_CUDA_MMQ_SRAM_LAYOUT_Q8_1,
     GGML_CUDA_MMQ_SRAM_LAYOUT_Q2_K,
     GGML_CUDA_MMQ_SRAM_LAYOUT_Q3_K,
@@ -133,6 +135,8 @@ static constexpr __host__ __device__ int ggml_cuda_mmq_get_sram_stride(ggml_cuda
     switch (sram_layout) {
         case GGML_CUDA_MMQ_SRAM_LAYOUT_Q8_0:
             return 2*MMQ_TILE_NE_K + 2*MMQ_TILE_NE_K/QI8_0 + 4;
+        case GGML_CUDA_MMQ_SRAM_LAYOUT_DT3:
+            return 4*MMQ_TILE_NE_K + 4*MMQ_TILE_NE_K/QI8_0 + 4;
         case GGML_CUDA_MMQ_SRAM_LAYOUT_Q8_1:
             return 2*MMQ_TILE_NE_K + 2*MMQ_TILE_NE_K/QI8_1 + 4;
         case GGML_CUDA_MMQ_SRAM_LAYOUT_Q2_K:
@@ -151,6 +155,7 @@ static constexpr __host__ __device__ int ggml_cuda_mmq_get_sram_stride(ggml_cuda
 }
 
 static_assert(ggml_cuda_mmq_get_sram_stride(GGML_CUDA_MMQ_SRAM_LAYOUT_Q8_0)  % 8 == 4, "Wrong padding.");
+static_assert(ggml_cuda_mmq_get_sram_stride(GGML_CUDA_MMQ_SRAM_LAYOUT_DT3)   % 8 == 4, "Wrong padding.");
 static_assert(ggml_cuda_mmq_get_sram_stride(GGML_CUDA_MMQ_SRAM_LAYOUT_Q8_1)  % 8 == 4, "Wrong padding.");
 static_assert(ggml_cuda_mmq_get_sram_stride(GGML_CUDA_MMQ_SRAM_LAYOUT_Q2_K)  % 8 == 4, "Wrong padding.");
 static_assert(ggml_cuda_mmq_get_sram_stride(GGML_CUDA_MMQ_SRAM_LAYOUT_Q3_K)  % 8 == 4, "Wrong padding.");
@@ -377,6 +382,7 @@ static constexpr __device__ int ggml_cuda_mmq_get_rows_per_warp(ggml_type type,
 #define MMQ_DP4A_TXS_Q8_0    tile_x_sizes{I*MMQ_TILE_NE_K*2 + I, I*MMQ_TILE_NE_K*2/QI8_0 + I/(QI8_0/2), 0}
 #define MMQ_DP4A_TXS_Q8_0_16 tile_x_sizes{I*MMQ_TILE_NE_K*2 + I, I*MMQ_TILE_NE_K*4/QI8_0 + I/(QI8_0/4), 0}
 #define MMQ_DP4A_TXS_Q8_1    tile_x_sizes{I*MMQ_TILE_NE_K*2 + I, I*MMQ_TILE_NE_K*2/QI8_1 + I/(QI8_1/2), 0}
+#define MMQ_DP4A_TXS_DT3     tile_x_sizes{I*MMQ_TILE_NE_K*4 + I, I*MMQ_TILE_NE_K*4/QI8_0 + I,           0}
 #define MMQ_DP4A_TXS_Q2_K    tile_x_sizes{I*MMQ_TILE_NE_K*2 + I, I*MMQ_TILE_NE_K         + I,           0}
 #define MMQ_DP4A_TXS_Q3_K    tile_x_sizes{I*MMQ_TILE_NE_K*2 + I, I,                                     I*MMQ_TILE_NE_K/8 + I/8}
 #define MMQ_DP4A_TXS_Q4_K    tile_x_sizes{I*MMQ_TILE_NE_K   + I, I*MMQ_TILE_NE_K/QI4_K,                 I*MMQ_TILE_NE_K/8 + I/8}
@@ -387,6 +393,7 @@ static constexpr __host__ __device__ tile_x_sizes mmq_get_dp4a_tile_x_sizes(ggml
     switch (type) {
         case GGML_TYPE_Q1_0:    return MMQ_DP4A_TXS_Q8_0;
         case GGML_TYPE_Q2_0:    return MMQ_DP4A_TXS_Q8_0;
+        case GGML_TYPE_DT3:     return MMQ_DP4A_TXS_DT3;
         case GGML_TYPE_Q4_0:    return MMQ_DP4A_TXS_Q4_0;
         case GGML_TYPE_Q4_1:    return MMQ_DP4A_TXS_Q4_1;
         case GGML_TYPE_Q5_0:    return MMQ_DP4A_TXS_Q8_0;
@@ -550,6 +557,12 @@ static constexpr __device__ ggml_cuda_mmq_util_funcs ggml_cuda_mmq_get_util_func
                     ggml_cuda_mmq_load_tiles_q2_0<type, J, fallback>,
                     ggml_cuda_mmq_vec_dot_q8_0_q8_1_dp4a<type, J, fallback>,
                     ggml_cuda_mmq_write_back_dp4a<type, J, fallback>);
+            case GGML_TYPE_DT3:
+                return ggml_cuda_mmq_util_funcs(
+                    VDR_Q8_0_Q8_1_MMQ,
+                    ggml_cuda_mmq_load_tiles_dt3<type, J, fallback>,
+                    ggml_cuda_mmq_vec_dot_dt3_q8_1_dp4a<type, J, fallback>,
+                    ggml_cuda_mmq_write_back_dp4a<type, J, fallback>);
             case GGML_TYPE_Q4_0:
                 return ggml_cuda_mmq_util_funcs(
                     VDR_Q4_0_Q8_1_MMQ,
@@ -714,6 +727,12 @@ static constexpr __device__ ggml_cuda_mmq_util_funcs ggml_cuda_mmq_get_util_func
                 ggml_cuda_mmq_load_tiles_q2_0<type, J, fallback>,
                 ggml_cuda_mmq_vec_dot_q8_0_q8_1_mma<type, J, fallback, MMQ_Q8_1_DS_LAYOUT_D4>,
                 ggml_cuda_mmq_write_back_mma<type, J, fallback>);
+        case GGML_TYPE_DT3:
+            return ggml_cuda_mmq_util_funcs(
+                -1,
+                ggml_cuda_mmq_load_tiles_dt3<type, J, fallback>,
+                ggml_cuda_mmq_vec_dot_dt3_q8_1_mma<type, J, fallback>,
+                ggml_cuda_mmq_write_back_mma<type, J, fallback>);
         case GGML_TYPE_Q4_0:
             return ggml_cuda_mmq_util_funcs(
                 -1,
@@ -1565,6 +1584,7 @@ void mul_mat_q_case(ggml_backend_cuda_context & ctx, const mmq_args & args, cuda
 
 extern DECL_MMQ_CASE(GGML_TYPE_Q1_0);
 extern DECL_MMQ_CASE(GGML_TYPE_Q2_0);
+extern DECL_MMQ_CASE(GGML_TYPE_DT3);
 extern DECL_MMQ_CASE(GGML_TYPE_Q4_0);
 extern DECL_MMQ_CASE(GGML_TYPE_Q4_1);
 extern DECL_MMQ_CASE(GGML_TYPE_Q5_0);

ggml/src/ggml-cuda/template-instances/generate_cu_files.py

@@ -37,6 +37,7 @@ SOURCE_FATTN_MMA_CASE = "DECL_FATTN_MMA_F16_CASE({head_size_kq}, {head_size_v},
 TYPES_MMQ = [
     "GGML_TYPE_Q1_0",
     "GGML_TYPE_Q2_0",
+    "GGML_TYPE_DT3",
     "GGML_TYPE_Q4_0", "GGML_TYPE_Q4_1", "GGML_TYPE_Q5_0", "GGML_TYPE_Q5_1", "GGML_TYPE_Q8_0",
     "GGML_TYPE_Q2_K", "GGML_TYPE_Q3_K", "GGML_TYPE_Q4_K", "GGML_TYPE_Q5_K", "GGML_TYPE_Q6_K",
     "GGML_TYPE_IQ2_XXS", "GGML_TYPE_IQ2_XS", "GGML_TYPE_IQ2_S", "GGML_TYPE_IQ3_XXS", "GGML_TYPE_IQ3_S",

ggml/src/ggml-cuda/template-instances/mmq-instance-dt3.cu

@@ -0,0 +1,5 @@
+// This file has been autogenerated by generate_cu_files.py, do not edit manually.
+
+#include "../mmq.cuh"
+
+DECL_MMQ_CASE(GGML_TYPE_DT3);

tests/test-dt3-gpu.cpp

@@ -18,7 +18,10 @@
 // the elementwise maximum is reported as information only, since it explodes
 // on cancellation whenever a true output value is near zero.
 //
-// MUL_MAT with more destination columns than the MMVQ limit falls back to
+// MUL_MAT with more destination columns than the MMVQ limit takes the MMQ
+// path where the backend implements it for DT3 (CUDA on Ampere-class
+// hardware and newer): integer dot products in the same numerical regime as
+// MMVQ, judged just as strictly. Backends without DT3 MMQ fall back to
 // dequantization + cuBLAS GEMM, which on fast-fp16 hardware rounds the
 // dequantized weights to fp16. DT3 weights (d1*t1 + d2*t2, the sum of two
 // fp16-scaled terms) are generally NOT fp16-representable, so that path is
@@ -264,7 +267,8 @@ static int test_mul_mat(ggml_backend_t backend, ggml_type type, const std::vecto
                         const std::vector<float> & ref_w, const std::vector<float> & y, bool strict) {
     int num_failed = 0;
 
-    const int ncols_dst[] = {1, 2, 5, 8, 16};
+    // 100 exercises a wide MMQ tile with a clamped last column block
+    const int ncols_dst[] = {1, 2, 5, 8, 16, 100};
 
     // the same weights as the fp16 GEMM fallback sees them
     std::vector<float> ref_w16(ref_w.size());
@@ -274,6 +278,8 @@ static int test_mul_mat(ggml_backend_t backend, ggml_type type, const std::vecto
 
     std::vector<std::vector<float>> results;
 
+    bool n16_integer = false; // whether the ncols_dst = 16 run took an integer (MMQ) path
+
     for (int c = 0; c < (int)(sizeof(ncols_dst)/sizeof(ncols_dst[0])); ++c) {
         const int n = ncols_dst[c];
 
@@ -324,25 +330,33 @@ static int test_mul_mat(ggml_backend_t backend, ggml_type type, const std::vecto
         const mat_err err16 = compare_mat(gpu, ref16);
 
         // n <= 8 is the MMVQ path with exact integer dot products, judged
-        // against the exact reference. Larger n is the dequantize + GEMM
-        // fallback whose numerics (fp16 or TF32 compute, depending on the
+        // against the exact reference (this mirrors MMVQ_MAX_BATCH_SIZE (8)
+        // from ggml-cuda/mmvq.cu by hand, because the constant and the
+        // per-arch should_use_mmvq tables are not exported). Larger n takes
+        // the MMQ path where the backend implements it for this type: integer
+        // dot products in the same numerical regime as MMVQ, judged just as
+        // strictly. Backends without MMQ for the type fall back to dequantize
+        // + GEMM, whose numerics (fp16 or TF32 compute, depending on the
         // hardware and on GGML_CUDA_CUBLAS_COMPUTE_TYPE) are cuBLAS's, not
-        // ours: for the strict type it is gated below by bit-identity with
-        // the same GEMM on an F16 tensor, and only reported here.
-        // This mirrors MMVQ_MAX_BATCH_SIZE (8) from ggml-cuda/mmvq.cu by hand,
-        // because the constant and the per-arch should_use_mmvq tables are not
-        // exported. If upstream raises the limit, or an architecture routes a
-        // larger batch through MMVQ, this gating goes stale silently: n = 16
-        // would take the MMVQ path but still be judged as the GEMM one, which
-        // only loosens the check, never tightens it. Whoever touches the MMVQ
-        // dispatch should revisit this line.
-        const bool   is_mmvq = n <= 8;
-        const bool   gated   = is_mmvq || !strict;
-        const double err_gate = is_mmvq ? err.norm_rel : (err.norm_rel < err16.norm_rel ? err.norm_rel : err16.norm_rel);
+        // ours: for the strict type that run is gated below by bit-identity
+        // with the same GEMM on an F16 tensor, and only reported here.
+        // The two regimes are told apart by the result itself: an integer path
+        // lands within float rounding of the exact reference, a fp16/TF32 GEMM
+        // stays orders of magnitude above it. A broken MMQ kernel cannot hide
+        // in the GEMM class: it would then have to be bit-identical to the F16
+        // GEMM control below, which an integer path never is.
+        const bool   is_mmvq  = n <= 8;
+        const bool   integer_path = is_mmvq || err.norm_rel <= 1e-5;
+        if (n == 16) {
+            n16_integer = integer_path;
+        }
+        const bool   gated   = integer_path || !strict;
+        const double err_gate = integer_path ? err.norm_rel : (err.norm_rel < err16.norm_rel ? err.norm_rel : err16.norm_rel);
         const double tol      = strict ? 1e-5 : 1e-2;
         const bool   failed   = gated && err_gate > tol;
-        printf("%s: %s mul_mat GPU, ncols_dst = %2d (%s): norm rel err vs exact ref = %g, vs fp16 ref = %g (max elem rel: %g)\n",
-               failed ? "FAILED" : gated ? "OK" : "INFO", ggml_type_name(type), n, is_mmvq ? "MMVQ" : "GEMM",
+        printf("%s: %s mul_mat GPU, ncols_dst = %3d (%s): norm rel err vs exact ref = %g, vs fp16 ref = %g (max elem rel: %g)\n",
+               failed ? "FAILED" : gated ? "OK" : "INFO", ggml_type_name(type), n,
+               is_mmvq ? "MMVQ" : integer_path ? "MMQ" : "GEMM",
                err.norm_rel, err16.norm_rel, err.max_rel);
         if (failed) {
             num_failed++;
@@ -352,24 +366,26 @@ static int test_mul_mat(ggml_backend_t backend, ggml_type type, const std::vecto
         ggml_free(ctx);
     }
 
-    // MMVQ vs the dequantization-based path: the first 8 columns of the GEMM
-    // run must match the ncols_dst = 8 MMVQ run to fp16 weight rounding
+    // MMVQ vs the batched path: the first 8 columns of the ncols_dst = 16 run
+    // must match the ncols_dst = 8 MMVQ run. When the batched run took the
+    // integer MMQ path both sides are exact to float rounding of the
+    // accumulation; against a GEMM fallback the gate is fp16 weight rounding.
     {
-        const std::vector<float> & mmvq = results[3]; // n = 8
-        const std::vector<float> & gemm = results[4]; // n = 16
+        const std::vector<float> & mmvq  = results[3]; // n = 8
+        const std::vector<float> & batch = results[4]; // n = 16
         double num = 0.0;
         double den = 0.0;
         for (int j = 0; j < 8; ++j) {
             for (int r = 0; r < NROWS; ++r) {
-                const double diff = (double)mmvq[(size_t)j*NROWS + r] - (double)gemm[(size_t)j*NROWS + r];
+                const double diff = (double)mmvq[(size_t)j*NROWS + r] - (double)batch[(size_t)j*NROWS + r];
                 num += diff*diff;
                 den += (double)mmvq[(size_t)j*NROWS + r]*(double)mmvq[(size_t)j*NROWS + r];
             }
         }
         const double norm_rel = sqrt(num/den);
-        const double tol      = strict ? 5e-3 : 1e-2;
-        printf("%s: %s MMVQ vs GEMM path on shared columns: norm rel err = %g\n",
-               norm_rel <= tol ? "OK" : "FAILED", ggml_type_name(type), norm_rel);
+        const double tol      = !strict ? 1e-2 : n16_integer ? 1e-5 : 5e-3;
+        printf("%s: %s MMVQ vs %s path on shared columns: norm rel err = %g\n",
+               norm_rel <= tol ? "OK" : "FAILED", ggml_type_name(type), n16_integer ? "MMQ" : "GEMM", norm_rel);
         if (norm_rel > tol) {
             num_failed++;
         }
@@ -383,8 +399,14 @@ static int test_mul_mat(ggml_backend_t backend, ggml_type type, const std::vecto
     // fallback at a different accumulator precision than its default F16
     // GEMM (Vulkan forces fp32 accumulators for DT3), so the F16 control is
     // run at both the default and the F32-forced precision and bit-identity
-    // with either one passes.
-    if (strict) {
+    // with either one passes. When the ncols_dst = 16 run took the integer
+    // MMQ path there is no dequantization involved and no GEMM to compare
+    // against — that run was already gated strictly above.
+    if (strict && n16_integer) {
+        printf("OK: %s ncols_dst = 16 took the integer MMQ path, F16 GEMM bit-identity control not applicable\n",
+               ggml_type_name(type));
+    }
+    if (strict && !n16_integer) {
         int    n_mismatch_best = -1;
         double max_diff_best   = 0.0;
 
@@ -522,7 +544,7 @@ int main(void) {
 
     // activations: integers with amax 127 in every 32-element chunk of every
     // column, so their q8_1 quantization is exact
-    std::vector<float> y((size_t)NCOLS*16);
+    std::vector<float> y((size_t)NCOLS*100);
     for (size_t i = 0; i < y.size(); ++i) {
         y[i] = i % 32 == 0 ? 127.0f : (float)((int)(rng_next() % 255) - 127);
     }