cuda : decode DT3 bytes once in the MMVQ vec_dot

The old vec_dot decoded every element with its own pair of multiplications
(256 inlined get_trit per block, each qs byte re-read 5 times). Decode each
byte once instead, iterating q -> (q*3) & 0xFF two bytes at a time in 16-bit
lanes, and accumulate dp4a over base-3 digits in {0, 1, 2}; one extra dp4a
with 0x01010101 per q8_1 int, shared by both planes, turns the digit sums
back into trit sums in exact integer arithmetic, so the result stays
bit-identical to the per-trit decode. The qh bytes keep their own 4-digit
path so the padding digit is never decoded.

ggml/src/ggml-cuda/vecdotq.cuh

@@ -776,26 +776,87 @@ static __device__ __forceinline__ float vec_dot_dt3_q8_1(
     // below are compile-time constants, so the decode folds into shifts and masks.
     GGML_UNUSED(iqs);
 
-    int sumi1[4] = {0, 0, 0, 0};
-    int sumi2[4] = {0, 0, 0, 0};
+    // The dot product is accumulated over base-3 digits in {0, 1, 2} instead of
+    // trits in {-1, 0, +1}: sum((digit - 1)*u) == sum(digit*u) - sum(u), and sum(u)
+    // is one extra dp4a with 0x01010101 shared by both planes. This allows decoding
+    // each packed byte once with the iteration q -> (q*3) & 0xFF, whose step n
+    // exposes base-3 digit n of the byte in bits 8..9 of q*3, instead of
+    // re-multiplying the byte by a power of 3 for every one of its 5 elements.
+    // Two bytes are iterated at a time in the 16-bit lanes of one int: a lane
+    // holds q < 256, so q*3 < 768 never carries into the neighbouring lane.
+    // The subtraction sumi - sumu happens in exact integer arithmetic, so the
+    // result is bit-identical to decoding the trits one by one.
+
+    int sumi[2][4] = {{0, 0, 0, 0}, {0, 0, 0, 0}}; // per plane, per q8_1 chunk: sum(digit*u)
+    int sumu[4]    = {0, 0, 0, 0};                 // per q8_1 chunk: sum(u)
 
 #pragma unroll
-    for (int j = 0; j < 4; ++j) {
+    for (int p = 0; p < 2; ++p) {
+        // qs[0..16): 4 quads of consecutive bytes, byte m holds elements m + 16*n
 #pragma unroll
-        for (int k = 0; k < 8; ++k) {
-            const int u = get_int_b4(bq8_1[j].qs, k);
+        for (int g = 0; g < 4; ++g) {
+            const int x = get_int_b4(bq_dt3->qs[p], g);
+            int qa = (x >> 0) & 0x00FF00FF; // bytes 4*g + 0 and 4*g + 2
+            int qb = (x >> 8) & 0x00FF00FF; // bytes 4*g + 1 and 4*g + 3
+#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;
 
-            int v1 = 0;
-            int v2 = 0;
-#pragma unroll
-            for (int l = 0; l < 4; ++l) {
-                const int i = 32*j + 4*k + l;
-                v1 |= (ggml_cuda_dt3_get_trit(bq_dt3->qs[0], bq_dt3->qh[0], i) & 0xFF) << (8*l);
-                v2 |= (ggml_cuda_dt3_get_trit(bq_dt3->qs[1], bq_dt3->qh[1], i) & 0xFF) << (8*l);
+                const int i = 16*n + 4*g; // first of the 4 consecutive elements
+                const int j = i / 32;
+                const int u = get_int_b4(bq8_1[j].qs, (i % 32)/4);
+                if (p == 0) {
+                    sumu[j] = ggml_cuda_dp4a(0x01010101, u, sumu[j]);
                 }
+                sumi[p][j] = ggml_cuda_dp4a(dig, u, sumi[p][j]);
+            }
+        }
+        // qs[16..24): 2 quads, byte 16 + m holds elements 80 + m + 8*n
+#pragma unroll
+        for (int g = 0; g < 2; ++g) {
+            const int x = get_int_b4(bq_dt3->qs[p], 4 + g);
+            int qa = (x >> 0) & 0x00FF00FF;
+            int qb = (x >> 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;
 
-            sumi1[j] = ggml_cuda_dp4a(v1, u, sumi1[j]);
-            sumi2[j] = ggml_cuda_dp4a(v2, u, sumi2[j]);
+                const int i = 80 + 8*n + 4*g;
+                const int j = i / 32;
+                const int u = get_int_b4(bq8_1[j].qs, (i % 32)/4);
+                if (p == 0) {
+                    sumu[j] = ggml_cuda_dp4a(0x01010101, u, sumu[j]);
+                }
+                sumi[p][j] = ggml_cuda_dp4a(dig, u, sumi[p][j]);
+            }
+        }
+        // qh[0..2): byte b holds elements 120 + b + 2*n for n = 0..4 — only 4
+        // digits are iterated, the 5th is packing padding and must never decode
+        {
+            int q = bq_dt3->qh[p][0] | (bq_dt3->qh[p][1] << 16);
+#pragma unroll
+            for (int s = 0; s < 2; ++s) { // one q8_1 int: elements 120 + 4*s .. 123 + 4*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;
+                }
+                const int u = get_int_b4(bq8_1[3].qs, 6 + s);
+                if (p == 0) {
+                    sumu[3] = ggml_cuda_dp4a(0x01010101, u, sumu[3]);
+                }
+                sumi[p][3] = ggml_cuda_dp4a(dig, u, sumi[p][3]);
+            }
         }
     }
 
@@ -806,7 +867,7 @@ static __device__ __forceinline__ float vec_dot_dt3_q8_1(
 #pragma unroll
     for (int j = 0; j < 4; ++j) {
         const float d8 = __low2float(bq8_1[j].ds);
-        sumf += d8 * (d1*sumi1[j] + d2*sumi2[j]);
+        sumf += d8 * (d1*(sumi[0][j] - sumu[j]) + d2*(sumi[1][j] - sumu[j]));
     }
     return sumf;
 }