diff --git a/server/CMakeLists.txt b/server/CMakeLists.txt index a969fa9d7..18eefd778 100644 --- a/server/CMakeLists.txt +++ b/server/CMakeLists.txt @@ -1143,6 +1143,33 @@ if(DFLASH27B_TESTS) endif() list(APPEND _raw_unit_test_targets test_rocmfp_mix_slice_matvec) endif() + if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/test/test_rocmfp_mix_mmq.cpp") + # Correctness gate for the BATCHED (MMQ) path of qtypes 105/106. Backend-general + # for the same reason as the tests below it, and doubly so here: the batched path + # was gated to RDNA behind an env var, so it had never been exercised on CUDA at + # all and had no correctness test on either vendor. + add_executable(test_rocmfp_mix_mmq test/test_rocmfp_mix_mmq.cpp) + if(DFLASH27B_GPU_BACKEND STREQUAL "hip") + set_source_files_properties(test/test_rocmfp_mix_mmq.cpp PROPERTIES LANGUAGE HIP) + set_target_properties(test_rocmfp_mix_mmq PROPERTIES HIP_ARCHITECTURES "${_dflash_archs}") + target_compile_definitions(test_rocmfp_mix_mmq PRIVATE GGML_USE_HIP) + else() + set_source_files_properties(test/test_rocmfp_mix_mmq.cpp PROPERTIES LANGUAGE CUDA) + set_target_properties(test_rocmfp_mix_mmq PROPERTIES CUDA_ARCHITECTURES "${_dflash_archs}") + endif() + target_include_directories(test_rocmfp_mix_mmq PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/deps/llama.cpp/ggml/include + ${CMAKE_CURRENT_SOURCE_DIR}/deps/llama.cpp/ggml/src/ggml-cuda) + target_link_libraries(test_rocmfp_mix_mmq PRIVATE + ggml ggml-base ${DFLASH27B_GGML_BACKEND_TARGET}) + if(DFLASH27B_GPU_BACKEND STREQUAL "cuda") + find_package(CUDAToolkit REQUIRED) + target_link_libraries(test_rocmfp_mix_mmq PRIVATE CUDA::cudart) + else() + target_link_libraries(test_rocmfp_mix_mmq PRIVATE hip::host) + endif() + list(APPEND _raw_unit_test_targets test_rocmfp_mix_mmq) + endif() if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/test/test_rocmfp_mix_gateup_glu.cpp") # Backend-general: the kernels declare cudaStream_t and ggml's vendors/hip.h maps the # cuda* spellings onto hip*, so the same source builds for both. These were hip-only, diff --git a/server/deps/llama.cpp/ggml/include/ggml-cuda.h b/server/deps/llama.cpp/ggml/include/ggml-cuda.h index ac2fc92d4..eb86c34b0 100644 --- a/server/deps/llama.cpp/ggml/include/ggml-cuda.h +++ b/server/deps/llama.cpp/ggml/include/ggml-cuda.h @@ -24,6 +24,20 @@ extern "C" { // wider batches remain on the MMQ path. #define GGML_CUDA_DS4_MIX_MMV_MAX_TOKENS 5 +// Batched (MMQ) path for the mix qtypes 105/106. Off unless a caller opts in: +// the measured benefit belongs to models whose mix tensors are dense and go +// through ggml_mul_mat, so the model backend that wants it enables it for +// itself rather than every model paying for one model's win. +// DFLASH_MIX_MMQ=0/1 (legacy: DFLASH_DS4_MIX_MMQ_PREFILL) overrides both. +GGML_BACKEND_API bool ggml_cuda_mix_mmq_enabled(void); +// Force the toggle regardless of environment, so one process can A/B both +// paths; test_rocmfp_mix_mmq uses this against the validated matvec kernel. +GGML_BACKEND_API void ggml_cuda_set_mix_mmq_enabled(bool enabled); +GGML_BACKEND_API void ggml_cuda_clear_mix_mmq_override(void); +// True when the environment pinned the toggle, so a backend's opt-in knows to +// leave an explicit operator choice alone. +GGML_BACKEND_API bool ggml_cuda_mix_mmq_env_pinned(void); + // backend API GGML_BACKEND_API ggml_backend_t ggml_backend_cuda_init(int device); diff --git a/server/deps/llama.cpp/ggml/src/ggml-cuda/mmq.cu b/server/deps/llama.cpp/ggml/src/ggml-cuda/mmq.cu index d782aaba1..cf2c96029 100644 --- a/server/deps/llama.cpp/ggml/src/ggml-cuda/mmq.cu +++ b/server/deps/llama.cpp/ggml/src/ggml-cuda/mmq.cu @@ -1,4 +1,5 @@ #include "common.cuh" +#include "ggml-cuda.h" #include "mmq.cuh" #include "quantize.cuh" #include "mmid.cuh" @@ -471,11 +472,130 @@ void ggml_cuda_op_mul_mat_q( GGML_UNUSED_VARS(src1, dst, src1_ddf_i, src1_padded_row_size); } +#ifndef GGML_CUDA_MIX_MMQ_DEFAULT +// OFF unless a caller opts in — ggml_cuda_set_mix_mmq_enabled() from a model +// backend, or DFLASH_MIX_MMQ=1. No in-tree backend opts in today, on purpose: +// this path applies to dense ggml_mul_mat, and DeepSeek-V4's 105/106 tensors +// are MoE experts reached through ggml_mul_mat_id, so turning it on changes +// nothing for DS4. That was MEASURED, not assumed — six serving configs +// byte-identical either way, a 3974-token prefill at 184.1 s +/-0.1%, and a +// controlled decode A/B at 21.72 / 21.72 / 21.74 tok/s on gfx1151. +// +// The win this path exists for was measured out of tree, on a dense 3.3 bpw +// consumer of these qtypes: 81.0 tok/s against 43.5 with it off, on the same +// 122-item task suite at equal quality. A process-wide default would extend +// that result to models it was never measured on, so the enable stays with +// whichever backend can show the measurement. +#define GGML_CUDA_MIX_MMQ_DEFAULT false +#endif + +// ── Batched path for the mix qtypes (105/106) ─────────────────────────── +// Without it, ne11 > 1 falls back to dequantize-to-bf16 + dense GEMM, which +// discards the whole point of a sub-4 bpw artifact for the duration of the +// multiply (measured: 48% of a 16-token speculative verify's GPU time inside +// dequantize_rocmfp{2,3}_mix_kernel). +// +// Runtime-settable rather than a function-local static so a single process can +// A/B both paths; test_rocmfp_mix_mmq compares them against the validated +// matvec kernel that way. Env var still provides the default. +static bool g_mix_mmq_forced = false; +static bool g_mix_mmq_forced_val = false; + +static const char * mix_mmq_env_value() { + const char * value = getenv("DFLASH_MIX_MMQ"); + if (value == nullptr) { + // Legacy spelling; the flag predates its use outside DS4 prefill. + value = getenv("DFLASH_DS4_MIX_MMQ_PREFILL"); + } + return value; +} + +bool ggml_cuda_mix_mmq_env_pinned() { + static const bool pinned = mix_mmq_env_value() != nullptr; + return pinned; +} + +// Precedence: an explicit set_ call (the test's A/B harness) beats the +// environment, which beats the compiled default. Backends opting in for their +// own model must therefore check ggml_cuda_mix_mmq_env_pinned() first, so an +// operator running DFLASH_MIX_MMQ=0 is not silently overridden. +bool ggml_cuda_mix_mmq_enabled() { + if (g_mix_mmq_forced) { + return g_mix_mmq_forced_val; + } + if (ggml_cuda_mix_mmq_env_pinned()) { + static const bool from_env = []() { + const char * value = mix_mmq_env_value(); + return !(value[0] == '0' && value[1] == '\0'); + }(); + return from_env; + } + return GGML_CUDA_MIX_MMQ_DEFAULT; +} + +void ggml_cuda_set_mix_mmq_enabled(bool enabled) { + g_mix_mmq_forced = true; + g_mix_mmq_forced_val = enabled; +} + +void ggml_cuda_clear_mix_mmq_override() { + g_mix_mmq_forced = false; +} + +// Widest ne11 at which the mix-qtype MMQ path still beats dequantize-to-bf16 + +// dense GEMM. MMQ is not uniformly better: it wins by a lot when the batch is +// narrow and loses when it is wide, because the dequant path hands a wide N to +// a well-tiled dense GEMM while the MMQ kernel does not tile as far. Measured +// out of tree on a dense 3.3 bpw artifact carrying 71 mix tensors (the only +// dense consumer of these qtypes so far), prefill tok/s, MMQ off -> on: +// +// ne11 8 16 64 256 1024 2048 +// gfx1151 1.80x 1.77x 1.65x 1.05x 0.89x 0.86x +// gfx1201 5.11x 4.02x 3.26x 1.83x 1.13x 0.98x +// +// so the crossover sits between 256 and 1024 on RDNA 3.5 and between 1024 and +// 2048 on RDNA 4. Decline MMQ past it rather than making the operator choose: +// decode (ne11 == 1) never reaches this gate at all, and a served request wants +// the narrow-batch win on its speculative-verify steps AND the wide-batch win +// on its prefill, which one process-wide boolean cannot deliver. +// +// NVIDIA has no width sweep yet — the H200 result behind this path (1.20x) was +// measured on ne11 4..16 verify batches. It therefore takes the conservative +// RDNA 3.5 bound, which keeps every measured NVIDIA win and declines only the +// widths nothing has measured there. Raise it with the env var once swept. +// +// CRITICAL: this cap is for the DENSE path only, and mul_mat_id opts out below. +// The cap is worth having because declining MMQ on a dense multiply hands the +// work to a well-tiled dense GEMM. Declining it on a mul_mat_id does something +// entirely different: it falls to the host-synchronised sort path (a +// cudaStreamSynchronize, a CPU-side id sort, an expand — and no CUDA-graph +// capture). Measured on H200 with 8 experts / top_k 2 at 512 tokens, +// q2_1_rocmfp2_mix: 0.412 ms with this cap applied vs 0.202 ms without, i.e. +// applying the dense cap to MoE costs 2.04x on exactly the widths a prefill +// uses. Same threshold, opposite consequence. +static int64_t mix_mmq_max_ne11(int cc) { + static const int64_t override_value = []() -> int64_t { + const char * value = getenv("DFLASH_MIX_MMQ_MAX_NE11"); + return value ? strtoll(value, nullptr, 10) : -1; + }(); + if (override_value >= 0) { + return override_value; + } + return GGML_CUDA_CC_IS_RDNA4(cc) ? 1024 : 256; +} + bool ggml_cuda_should_use_mmq(enum ggml_type type, int cc, int64_t ne11, int64_t n_experts) { #ifdef GGML_CUDA_FORCE_CUBLAS return false; #endif // GGML_CUDA_FORCE_CUBLAS + // Which dispatch is asking? The dense call sites pass 0 (or 1 for a 2-D + // src0 through the fused gate/up pair); ggml_cuda_mul_mat_id passes the + // real expert count. Only the mix qtypes act on this, and only to skip a + // width cap that would send MoE to the host-sync fallback — see + // mix_mmq_max_ne11. + const bool is_mul_mat_id = n_experts > 1; + bool mmq_supported; switch (type) { @@ -517,21 +637,49 @@ bool ggml_cuda_should_use_mmq(enum ggml_type type, int cc, int64_t ne11, int64_t GGML_CUDA_CC_IS_RDNA4(cc); break; case GGML_TYPE_Q2_1_ROCMFP2_MIX: { - static const bool mix_mmq_enabled = []() { - const char * value = getenv("DFLASH_DS4_MIX_MMQ_PREFILL"); - return value != nullptr && !(value[0] == '0' && value[1] == '\0'); - }(); - mmq_supported = mix_mmq_enabled && - (GGML_CUDA_CC_IS_RDNA3_5(cc) || GGML_CUDA_CC_IS_RDNA4(cc)); + // Batched path for the mix qtypes. Without it, ne11 > 1 falls back + // to dequantize-to-bf16 + dense GEMM, which throws away the whole + // point of a 3.3 bpw artifact for the duration of the multiply: + // measured 48% of a 16-token speculative verify's GPU time inside + // dequantize_rocmfp{2,3}_mix_kernel. + // + // Enabling it on gfx1201 measured 1.7-1.9x on that verify AND + // halved the batch-vs-sequential logit drift (0.774 -> 0.380), + // because the dequant path rounds through bf16 where MMQ keeps + // integer dot products. Faster and more faithful. + // + // NVIDIA is opted in behind the same env var so the claim can be + // measured rather than assumed: the DP4A tile these types declare + // is portable, but nothing has gated it here yet. + // + // Width-gated: see mix_mmq_max_ne11 for why wide batches decline. + mmq_supported = ggml_cuda_mix_mmq_enabled() && + (is_mul_mat_id || ne11 <= mix_mmq_max_ne11(cc)) && + (GGML_CUDA_CC_IS_RDNA3_5(cc) || GGML_CUDA_CC_IS_RDNA4(cc) || + GGML_CUDA_CC_IS_NVIDIA(cc)); break; } case GGML_TYPE_Q3_1_ROCMFP3_MIX: { - static const bool mix_mmq_enabled = []() { - const char * value = getenv("DFLASH_DS4_MIX_MMQ_PREFILL"); - return value != nullptr && !(value[0] == '0' && value[1] == '\0'); - }(); - mmq_supported = mix_mmq_enabled && - (GGML_CUDA_CC_IS_RDNA3_5(cc) || GGML_CUDA_CC_IS_RDNA4(cc)); + // Batched path for the mix qtypes. Without it, ne11 > 1 falls back + // to dequantize-to-bf16 + dense GEMM, which throws away the whole + // point of a 3.3 bpw artifact for the duration of the multiply: + // measured 48% of a 16-token speculative verify's GPU time inside + // dequantize_rocmfp{2,3}_mix_kernel. + // + // Enabling it on gfx1201 measured 1.7-1.9x on that verify AND + // halved the batch-vs-sequential logit drift (0.774 -> 0.380), + // because the dequant path rounds through bf16 where MMQ keeps + // integer dot products. Faster and more faithful. + // + // NVIDIA is opted in behind the same env var so the claim can be + // measured rather than assumed: the DP4A tile these types declare + // is portable, but nothing has gated it here yet. + // + // Width-gated: see mix_mmq_max_ne11 for why wide batches decline. + mmq_supported = ggml_cuda_mix_mmq_enabled() && + (is_mul_mat_id || ne11 <= mix_mmq_max_ne11(cc)) && + (GGML_CUDA_CC_IS_RDNA3_5(cc) || GGML_CUDA_CC_IS_RDNA4(cc) || + GGML_CUDA_CC_IS_NVIDIA(cc)); break; } default: diff --git a/server/deps/llama.cpp/ggml/src/ggml-cuda/rocmfp2_mix.cu b/server/deps/llama.cpp/ggml/src/ggml-cuda/rocmfp2_mix.cu index 342e396fc..045d0aba7 100644 --- a/server/deps/llama.cpp/ggml/src/ggml-cuda/rocmfp2_mix.cu +++ b/server/deps/llama.cpp/ggml/src/ggml-cuda/rocmfp2_mix.cu @@ -327,6 +327,40 @@ __device__ __forceinline__ float mix_warp_shfl_down(float v, int off) { #endif } +// Stage a block's 32 activations into registers with float4 loads. +// +// Measured motivation (H200, dense decode of these qtypes): the per-block j loop +// below reads xc[col0 + 0 .. 31], and nvcc emitted one 32-bit LDG per element — +// the overwhelming majority of the kernel's global load instructions, against a +// handful for the (already wide-staged) weights. ncu put the sibling qtype-105 +// kernel at 1.65% DRAM throughput, 97% L1 throughput and 98.4% L1 hit rate with +// 56% of warp cycles stalled on LG throttle: not bandwidth bound, bound purely +// on the NUMBER of narrow global load instructions. The 32 floats are contiguous +// (128 B), so four at a time costs one instruction instead of four. +// +// Bit-exact: this changes only how the SAME 32 f32 values reach registers. The +// j loop still consumes them in ascending order into the same acc chain. +// +// The alignment test is warp- AND block-uniform (col0 is always a multiple of +// MIX_QK, i.e. 128 B, so xc+col0 has xc's alignment), so the branch costs no +// divergence. The scalar arm is not dead code: ggml hands this kernel a src1 +// row pointer, and a strided/offset column would land unaligned. +__device__ __forceinline__ void mix_load_x32( + const float * __restrict__ xp, float (&xv)[MIX_QK]) { + if ((((uintptr_t) xp) & 15u) == 0) { + const float4 * __restrict__ p4 = (const float4 *) xp; + #pragma unroll + for (int k = 0; k < MIX_QK / 4; ++k) { + const float4 v = p4[k]; + xv[4*k + 0] = v.x; xv[4*k + 1] = v.y; + xv[4*k + 2] = v.z; xv[4*k + 3] = v.w; + } + } else { + #pragma unroll + for (int j = 0; j < MIX_QK; ++j) xv[j] = xp[j]; + } +} + // Accumulate one block's 32 terms directly into acc, in fixed j order, exactly // as the un-refactored loop did (acc += s * w_j * x_j). Adding each term into // the shared running acc — rather than forming a per-block partial sum first — @@ -335,8 +369,8 @@ __device__ __forceinline__ float mix_warp_shfl_down(float v, int off) { // tree reduction changes the rounding and flips tokens.) The block's byte loads // do not depend on acc, so unrolling the caller over several blocks lets the // compiler overlap their loads even though the acc-add chain stays serial. -__device__ __forceinline__ void mix_block_accum( - const uint8_t * __restrict__ b, const float * __restrict__ xc, int col0, +__device__ __forceinline__ void mix_block_accum_x( + const uint8_t * __restrict__ b, const float (&xv)[MIX_QK], int mode, const float * __restrict__ lut, float & acc) { // Stage the whole 10-byte block into registers with one 2-byte-wide copy, // then decode the fp2 codes out of registers instead of re-reading the @@ -386,7 +420,7 @@ __device__ __forceinline__ void mix_block_accum( #pragma unroll for (int j = 0; j < MIX_QK; ++j) { const float s = (j < MIX_QK/2) ? s0 : s1; - acc += s * mix_fp2_fixed(mix_fp2_code_u64(codes, j)) * xc[col0 + j]; + acc += s * mix_fp2_fixed(mix_fp2_code_u64(codes, j)) * xv[j]; } } else { // `lut` is the expert's 2 * MIX_K codebook ALREADY widened to f32 and staged in @@ -407,11 +441,41 @@ __device__ __forceinline__ void mix_block_accum( for (int j = 0; j < MIX_QK; ++j) { const float s = (j < MIX_QK/2) ? s0 : s1; const float * bk = (j < MIX_QK/2) ? bk0 : bk1; - acc += s * bk[mix_fp2_code_u64(codes, j)] * xc[col0 + j]; + acc += s * bk[mix_fp2_code_u64(codes, j)] * xv[j]; } } } +// Original (col0-addressed) entry point, kept for the MoE and 3-D slice kernels. +// Stages the activations itself; identical arithmetic. +__device__ __forceinline__ void mix_block_accum( + const uint8_t * __restrict__ b, const float * __restrict__ xc, int col0, + int mode, const float * __restrict__ lut, float & acc) { + float xv[MIX_QK]; + mix_load_x32(xc + col0, xv); + mix_block_accum_x(b, xv, mode, lut, acc); +} + +// Two output rows, ONE activation stage. The 32 activations a block consumes are +// the same for every row, so loading them once and folding both rows against the +// registers halves the (dominant) activation load issue per unit work. This is +// the same register-blocking the MoE kernel below already does — the dense 2-D +// kernel never got it because the model this family was built for is MoE and +// never took this path. +// +// Bit-exact: each row keeps its own accumulator and its own ascending-j fold over +// the same block sequence, so acc0/acc1 match the one-row-per-warp kernel element +// for element. +__device__ __forceinline__ void mix_block_accum2( + const uint8_t * __restrict__ b0, const uint8_t * __restrict__ b1, + const float * __restrict__ xp, int mode, const float * __restrict__ lut, + float & acc0, float & acc1) { + float xv[MIX_QK]; + mix_load_x32(xp, xv); + mix_block_accum_x(b0, xv, mode, lut, acc0); + mix_block_accum_x(b1, xv, mode, lut, acc1); +} + // The lane's block loop is unrolled by MIX_UNROLL into a SINGLE accumulator kept // in the exact original block order (acc += dot(blk), stride MIX_WARP), so the // f32 output is bit-for-bit identical to the un-unrolled path — required because @@ -428,7 +492,8 @@ __global__ void mix_matvec_rocmfp2_kernel( float * __restrict__ y, int in, int out, int64_t x_col_stride, int64_t y_col_stride) { const int warps_per_block = blockDim.x / MIX_WARP; - const int row = blockIdx.x * warps_per_block + (threadIdx.x / MIX_WARP); + const int warp = blockIdx.x * warps_per_block + (threadIdx.x / MIX_WARP); + const int row = warp * 2; // two output rows per warp const int lane = threadIdx.x % MIX_WARP; const int col = blockIdx.y; // Widen the 2 * MIX_K bf16 codebook to f32 in LDS once per workgroup instead of @@ -443,8 +508,14 @@ __global__ void mix_matvec_rocmfp2_kernel( if (row >= out) return; const int mode = (int) mode_ptr[0]; const int nb = in / MIX_QK; - const uint8_t * rowbase = data + (int64_t) row * nb * MIX_BLOCK_BYTES; - const float * xc = x + (int64_t) col * x_col_stride; + // `two` is warp-uniform, so the hot loop never diverges. It is false only for + // the tail warp of an odd `out`; that warp reuses row0's base so every load + // stays in bounds and acc1 is simply never stored. + const bool two = (row + 1) < out; + const uint8_t * rowbase0 = data + (int64_t) row * nb * MIX_BLOCK_BYTES; + const uint8_t * rowbase1 = two ? data + (int64_t) (row + 1) * nb * MIX_BLOCK_BYTES + : rowbase0; + const float * xc = x + (int64_t) col * x_col_stride; // f32 accumulate of f32-dequantized weights * f32 activations. The dequant // is bit-exact vs the reference (validated in ~/p4-validate/hip). This is // slightly higher precision than the f16 dequant->cuBLAS fallback; on the @@ -452,7 +523,7 @@ __global__ void mix_matvec_rocmfp2_kernel( // noise band, with every divergence a shared model-limited miss, a harness // answer-extraction artifact, or an HE formatting coin-flip -- no reasoning // regression. Kept f32 for simplicity/speed (no per-weight rounding ops). - float acc = 0.0f; + float acc0 = 0.0f, acc1 = 0.0f; int blk = lane; // Main body: MIX_UNROLL blocks per iteration, each accumulated in stride // order into the single acc. The blocks' byte loads are independent (only @@ -462,16 +533,27 @@ __global__ void mix_matvec_rocmfp2_kernel( #pragma unroll for (int u = 0; u < MIX_UNROLL; ++u) { const int b = blk + u * MIX_WARP; - mix_block_accum(rowbase + (int64_t) b * MIX_BLOCK_BYTES, xc, b * MIX_QK, mode, s_lut, acc); + mix_block_accum2(rowbase0 + (int64_t) b * MIX_BLOCK_BYTES, + rowbase1 + (int64_t) b * MIX_BLOCK_BYTES, + xc + b * MIX_QK, mode, s_lut, acc0, acc1); } } // Remainder: fewer than MIX_UNROLL strided blocks left for this lane. for (; blk < nb; blk += MIX_WARP) { - mix_block_accum(rowbase + (int64_t) blk * MIX_BLOCK_BYTES, xc, blk * MIX_QK, mode, s_lut, acc); + mix_block_accum2(rowbase0 + (int64_t) blk * MIX_BLOCK_BYTES, + rowbase1 + (int64_t) blk * MIX_BLOCK_BYTES, + xc + blk * MIX_QK, mode, s_lut, acc0, acc1); } #pragma unroll - for (int off = MIX_WARP/2; off > 0; off >>= 1) acc += mix_warp_shfl_down(acc, off); - if (lane == 0) y[(int64_t) col * y_col_stride + row] = acc; + for (int off = MIX_WARP/2; off > 0; off >>= 1) { + acc0 += mix_warp_shfl_down(acc0, off); + acc1 += mix_warp_shfl_down(acc1, off); + } + if (lane == 0) { + const int64_t o = (int64_t) col * y_col_stride + row; + y[o] = acc0; + if (two) y[o + 1] = acc1; + } } // ---- stream-sync-free fused MoE matvec (mul_mat_id) ---- @@ -878,7 +960,9 @@ bool ggml_cuda_rocmfp2_mix_mul_mat_vec( // give the scheduler finer packing/tail balance on this BW-bound matvec. const int warps_per_block = 2; // 64 threads const int threads = warps_per_block * MIX_WARP; - dim3 grid((out + warps_per_block - 1) / warps_per_block, ncols, 1); + // Each warp now owns TWO rows, so a block covers 2 * warps_per_block of them. + const int rows_per_block = 2 * warps_per_block; + dim3 grid((out + rows_per_block - 1) / rows_per_block, ncols, 1); mix_matvec_rocmfp2_kernel<<>>( (const uint8_t *) vx, book, mode_ptr, x, y, in, out, x_col_stride, y_col_stride); return true; diff --git a/server/deps/llama.cpp/ggml/src/ggml-cuda/rocmfp3_mix.cu b/server/deps/llama.cpp/ggml/src/ggml-cuda/rocmfp3_mix.cu index f24031509..017bb405e 100644 --- a/server/deps/llama.cpp/ggml/src/ggml-cuda/rocmfp3_mix.cu +++ b/server/deps/llama.cpp/ggml/src/ggml-cuda/rocmfp3_mix.cu @@ -297,8 +297,42 @@ __device__ __forceinline__ float mix_warp_shfl_down(float v, int off) { // tree reduction changes the rounding and flips tokens.) The block's byte loads // do not depend on acc, so unrolling the caller over several blocks lets the // compiler overlap their loads even though the acc-add chain stays serial. -__device__ __forceinline__ void mix_block_accum( - const uint8_t * __restrict__ b, const float * __restrict__ xc, int col0, +// Stage a block's 32 activations into registers with float4 loads. +// +// Measured motivation (H200, dense decode of these qtypes): the per-block j loop +// below reads xc[col0 + 0 .. 31], and nvcc emitted 32 separate 32-bit +// LDG.E.CONSTANT for them — 320 of the kernel's 355 global load instructions, +// against 35 for the weights. ncu put the kernel at 1.65% DRAM throughput, +// 97% L1 throughput and 98.4% L1 hit rate with 56% of all warp cycles stalled +// on LG throttle: not bandwidth bound, bound purely on the NUMBER of narrow +// global load instructions. The 32 floats are contiguous (128 B), so four at a +// time costs one instruction instead of four. +// +// Bit-exact: this changes only how the SAME 32 f32 values reach registers. The +// j loop still consumes them in ascending order into the same acc chain. +// +// The alignment test is warp- AND block-uniform (col0 is always a multiple of +// MIX_QK, i.e. 128 B, so xc+col0 has xc's alignment), so the branch costs no +// divergence. The scalar arm is not dead code: ggml hands this kernel a src1 +// row pointer, and a strided/offset column would land unaligned. +__device__ __forceinline__ void mix_load_x32( + const float * __restrict__ xp, float (&xv)[MIX_QK]) { + if ((((uintptr_t) xp) & 15u) == 0) { + const float4 * __restrict__ p4 = (const float4 *) xp; + #pragma unroll + for (int k = 0; k < MIX_QK / 4; ++k) { + const float4 v = p4[k]; + xv[4*k + 0] = v.x; xv[4*k + 1] = v.y; + xv[4*k + 2] = v.z; xv[4*k + 3] = v.w; + } + } else { + #pragma unroll + for (int j = 0; j < MIX_QK; ++j) xv[j] = xp[j]; + } +} + +__device__ __forceinline__ void mix_block_accum_x( + const uint8_t * __restrict__ b, const float (&xv)[MIX_QK], int mode, const float * __restrict__ lut, float & acc) { // Stage the whole 14-byte block into registers with one 2-byte-wide copy, // then decode the fp3 codes out of registers instead of re-reading the @@ -319,7 +353,7 @@ __device__ __forceinline__ void mix_block_accum( #pragma unroll for (int j = 0; j < MIX_QK; ++j) { const float s = (j < MIX_QK/2) ? s0 : s1; - acc += s * mix_fp3_fixed(mix_fp3_code(buf, j)) * xc[col0 + j]; + acc += s * mix_fp3_fixed(mix_fp3_code(buf, j)) * xv[j]; } } else { const float s0 = mix_ue4m3(m0 & 0x7F), s1 = mix_ue4m3(m1 & 0x7F); @@ -336,11 +370,41 @@ __device__ __forceinline__ void mix_block_accum( for (int j = 0; j < MIX_QK; ++j) { const float s = (j < MIX_QK/2) ? s0 : s1; const float * bk = (j < MIX_QK/2) ? bk0 : bk1; - acc += s * bk[mix_fp3_code(buf, j)] * xc[col0 + j]; + acc += s * bk[mix_fp3_code(buf, j)] * xv[j]; } } } +// Original (col0-addressed) entry point, kept for the MoE and 3-D slice kernels. +// Stages the activations itself; identical arithmetic. +__device__ __forceinline__ void mix_block_accum( + const uint8_t * __restrict__ b, const float * __restrict__ xc, int col0, + int mode, const float * __restrict__ lut, float & acc) { + float xv[MIX_QK]; + mix_load_x32(xc + col0, xv); + mix_block_accum_x(b, xv, mode, lut, acc); +} + +// Two output rows, ONE activation stage. The 32 activations a block consumes are +// the same for every row, so loading them once and folding both rows against the +// registers halves the (dominant) activation load issue per unit work. This is +// the same register-blocking the MoE kernel below already does — the dense 2-D +// kernel never got it because the model this family was built for is MoE and +// never took this path. +// +// Bit-exact: each row keeps its own accumulator and its own ascending-j fold over +// the same block sequence, so acc0/acc1 match the one-row-per-warp kernel element +// for element. +__device__ __forceinline__ void mix_block_accum2( + const uint8_t * __restrict__ b0, const uint8_t * __restrict__ b1, + const float * __restrict__ xp, int mode, const float * __restrict__ lut, + float & acc0, float & acc1) { + float xv[MIX_QK]; + mix_load_x32(xp, xv); + mix_block_accum_x(b0, xv, mode, lut, acc0); + mix_block_accum_x(b1, xv, mode, lut, acc1); +} + // The lane's block loop is unrolled by MIX_UNROLL into a SINGLE accumulator kept // in the exact original block order (acc += dot(blk), stride MIX_WARP), so the // f32 output is bit-for-bit identical to the un-unrolled path — required because @@ -357,7 +421,8 @@ __global__ void mix_matvec_rocmfp3_kernel( float * __restrict__ y, int in, int out, int64_t x_col_stride, int64_t y_col_stride) { const int warps_per_block = blockDim.x / MIX_WARP; - const int row = blockIdx.x * warps_per_block + (threadIdx.x / MIX_WARP); + const int warp = blockIdx.x * warps_per_block + (threadIdx.x / MIX_WARP); + const int row = warp * 2; // two output rows per warp const int lane = threadIdx.x % MIX_WARP; const int col = blockIdx.y; // Widen the 2 * MIX_K bf16 codebook to f32 in LDS once per workgroup rather than @@ -372,8 +437,14 @@ __global__ void mix_matvec_rocmfp3_kernel( if (row >= out) return; const int mode = (int) mode_ptr[0]; const int nb = in / MIX_QK; - const uint8_t * rowbase = data + (int64_t) row * nb * MIX_BLOCK_BYTES; - const float * xc = x + (int64_t) col * x_col_stride; + // `two` is warp-uniform, so the hot loop never diverges. It is false only for + // the tail warp of an odd `out`; that warp reuses row0's base so every load + // stays in bounds and acc1 is simply never stored. + const bool two = (row + 1) < out; + const uint8_t * rowbase0 = data + (int64_t) row * nb * MIX_BLOCK_BYTES; + const uint8_t * rowbase1 = two ? data + (int64_t) (row + 1) * nb * MIX_BLOCK_BYTES + : rowbase0; + const float * xc = x + (int64_t) col * x_col_stride; // f32 accumulate of f32-dequantized weights * f32 activations. The dequant // is bit-exact vs the reference (validated in ~/p4-validate/hip). This is // slightly higher precision than the f16 dequant->cuBLAS fallback; on the @@ -381,7 +452,7 @@ __global__ void mix_matvec_rocmfp3_kernel( // noise band, with every divergence a shared model-limited miss, a harness // answer-extraction artifact, or an HE formatting coin-flip -- no reasoning // regression. Kept f32 for simplicity/speed (no per-weight rounding ops). - float acc = 0.0f; + float acc0 = 0.0f, acc1 = 0.0f; int blk = lane; // Main body: MIX_UNROLL blocks per iteration, each accumulated in stride // order into the single acc. The blocks' byte loads are independent (only @@ -390,18 +461,35 @@ __global__ void mix_matvec_rocmfp3_kernel( for (; blk + 3 * MIX_WARP < nb; blk += MIX_UNROLL * MIX_WARP) { const int b0 = blk, b1 = blk + MIX_WARP; const int b2 = blk + 2 * MIX_WARP, b3 = blk + 3 * MIX_WARP; - mix_block_accum(rowbase + (int64_t) b0 * MIX_BLOCK_BYTES, xc, b0 * MIX_QK, mode, s_lut, acc); - mix_block_accum(rowbase + (int64_t) b1 * MIX_BLOCK_BYTES, xc, b1 * MIX_QK, mode, s_lut, acc); - mix_block_accum(rowbase + (int64_t) b2 * MIX_BLOCK_BYTES, xc, b2 * MIX_QK, mode, s_lut, acc); - mix_block_accum(rowbase + (int64_t) b3 * MIX_BLOCK_BYTES, xc, b3 * MIX_QK, mode, s_lut, acc); + mix_block_accum2(rowbase0 + (int64_t) b0 * MIX_BLOCK_BYTES, + rowbase1 + (int64_t) b0 * MIX_BLOCK_BYTES, + xc + b0 * MIX_QK, mode, s_lut, acc0, acc1); + mix_block_accum2(rowbase0 + (int64_t) b1 * MIX_BLOCK_BYTES, + rowbase1 + (int64_t) b1 * MIX_BLOCK_BYTES, + xc + b1 * MIX_QK, mode, s_lut, acc0, acc1); + mix_block_accum2(rowbase0 + (int64_t) b2 * MIX_BLOCK_BYTES, + rowbase1 + (int64_t) b2 * MIX_BLOCK_BYTES, + xc + b2 * MIX_QK, mode, s_lut, acc0, acc1); + mix_block_accum2(rowbase0 + (int64_t) b3 * MIX_BLOCK_BYTES, + rowbase1 + (int64_t) b3 * MIX_BLOCK_BYTES, + xc + b3 * MIX_QK, mode, s_lut, acc0, acc1); } // Remainder: fewer than MIX_UNROLL strided blocks left for this lane. for (; blk < nb; blk += MIX_WARP) { - mix_block_accum(rowbase + (int64_t) blk * MIX_BLOCK_BYTES, xc, blk * MIX_QK, mode, s_lut, acc); + mix_block_accum2(rowbase0 + (int64_t) blk * MIX_BLOCK_BYTES, + rowbase1 + (int64_t) blk * MIX_BLOCK_BYTES, + xc + blk * MIX_QK, mode, s_lut, acc0, acc1); } #pragma unroll - for (int off = MIX_WARP/2; off > 0; off >>= 1) acc += mix_warp_shfl_down(acc, off); - if (lane == 0) y[(int64_t) col * y_col_stride + row] = acc; + for (int off = MIX_WARP/2; off > 0; off >>= 1) { + acc0 += mix_warp_shfl_down(acc0, off); + acc1 += mix_warp_shfl_down(acc1, off); + } + if (lane == 0) { + const int64_t o = (int64_t) col * y_col_stride + row; + y[o] = acc0; + if (two) y[o + 1] = acc1; + } } // ---- stream-sync-free fused MoE matvec (mul_mat_id) ---- @@ -786,7 +874,9 @@ bool ggml_cuda_rocmfp3_mix_mul_mat_vec( // give the scheduler finer packing/tail balance on this BW-bound matvec. const int warps_per_block = 2; // 64 threads const int threads = warps_per_block * MIX_WARP; - dim3 grid((out + warps_per_block - 1) / warps_per_block, ncols, 1); + // Each warp now owns TWO rows, so a block covers 2 * warps_per_block of them. + const int rows_per_block = 2 * warps_per_block; + dim3 grid((out + rows_per_block - 1) / rows_per_block, ncols, 1); mix_matvec_rocmfp3_kernel<<>>( (const uint8_t *) vx, book, mode_ptr, x, y, in, out, x_col_stride, y_col_stride); return true; diff --git a/server/test/test_rocmfp_mix_mmq.cpp b/server/test/test_rocmfp_mix_mmq.cpp new file mode 100644 index 000000000..a44016180 --- /dev/null +++ b/server/test/test_rocmfp_mix_mmq.cpp @@ -0,0 +1,522 @@ +// Correctness gate for the BATCHED (MMQ) path of the mix qtypes 105/106. +// +// Why this exists. With `ne11 > 1` these qtypes had no reachable batched +// kernel: `ggml_cuda_should_use_mmq` gated MMQ to RDNA behind an env var, so +// every batched multiply fell back to dequantize-to-bf16 + dense GEMM. That is +// correct but throws the format away for the duration of the multiply +// (measured: 48% of a 16-token speculative verify's GPU time inside +// dequantize_rocmfp{2,3}_mix_kernel). Turning MMQ on is worth 1.7-1.9x on +// gfx1201 — but a default cannot flip on a throughput number alone, and there +// was NO correctness test for mix MMQ anywhere. This is that test. +// +// What is compared, and why not bit-identity. MMQ quantizes the ACTIVATIONS to +// q8_1 before the dot product; the dequant path keeps them in f32 and rounds +// the WEIGHTS through bf16 instead. The two are different algorithms, so they +// cannot be bit-identical and demanding that would be wrong. What can be +// demanded is that MMQ is **no worse than the path it replaces**, both measured +// against the same reference: +// +// reference = ggml_cuda_rocmfp{2,3}_mix_mul_mat_vec, the already-validated +// kernel that batch-1 decode uses and that the greedy-output +// correctness gate hashes. +// +// So the assertion is err(MMQ) <= tolerance * err(dequant). If MMQ were +// decoding against the wrong codebook, mis-striding the tile, or ignoring the +// per-tensor mode byte, its error would be O(magnitude) and this fails loudly, +// which is the failure mode that matters: wrong numbers, not a crash. +// +// Coverage is deliberately across the axes that differ between the two callers +// of these qtypes: +// - both qtypes (105 = 3-bit mix, 106 = 2-bit mix) +// - both codebook modes (0 = fixed levels, 1 = learned codebook) +// - several `ne11` including 1 (matvec territory) and 16 (the speculative +// verify width) and 64 (prefill territory) +// - n_experts 1 (dense mul_mat) and > 1 (DS4, MoE) — a kernel that ignored the +// expert stride would pass the dense case and fail here + +#include "ds4_test_gpu_runtime.h" + +#include "ggml.h" +#include "ggml-alloc.h" +#include "ggml-backend.h" +#include "ggml-cuda.h" + +#include +#include +#include +#include +#include +#include +#include + +// ggml-cuda.h also declares the mix-MMQ runtime override, so one process can +// A/B both paths against the same weights; without it the env var is read once +// per process and the comparison would need two runs. + +extern "C" void ggml_cuda_rocmfp3_mix_register_host( + const void * base, size_t nb02, int n_experts, int out, int in, + const void * codebooks_bf16_host, const uint8_t * modes_host, + const uint8_t * rotations_host); +extern "C" void ggml_cuda_rocmfp3_mix_unregister(const void * base); +extern "C" void ggml_cuda_rocmfp2_mix_register_host( + const void * base, size_t nb02, int n_experts, int out, int in, + const void * codebooks_bf16_host, const uint8_t * modes_host, + const uint8_t * rotations_host); +extern "C" void ggml_cuda_rocmfp2_mix_unregister(const void * base); + +bool ggml_cuda_rocmfp3_mix_mul_mat_vec( + const void * vx, const float * x, float * y, int in, int out, int ncols, + int64_t x_col_stride, int64_t y_col_stride, cudaStream_t stream); +bool ggml_cuda_rocmfp2_mix_mul_mat_vec( + const void * vx, const float * x, float * y, int in, int out, int ncols, + int64_t x_col_stride, int64_t y_col_stride, cudaStream_t stream); + +namespace { + +int g_fails = 0; + +void fail(const std::string & msg) { + std::fprintf(stderr, "FAIL: %s\n", msg.c_str()); + ++g_fails; +} + +uint32_t g_rng = 0x2545F491u; +uint32_t rnd() { + g_rng ^= g_rng << 13; g_rng ^= g_rng >> 17; g_rng ^= g_rng << 5; + return g_rng; +} + +struct MixType { + ggml_type type; + const char * label; + void (*register_host)(const void *, size_t, int, int, int, const void *, + const uint8_t *, const uint8_t *); + void (*unregister)(const void *); + bool (*matvec)(const void *, const float *, float *, int, int, int, + int64_t, int64_t, cudaStream_t); +}; + +// Fill plausible blocks. The encoder is not reimplemented here: the kernels +// decode whatever bytes are present and this test compares kernel against +// kernel, so arbitrary-but-valid bytes exercise the same decode paths. +// Scale indices are kept inside the finite UE4M3 range — >0x7E decodes to 0.0 +// and whole half-blocks would vanish, weakening the comparison. +std::vector make_blocks(size_t bytes, size_t block_bytes) { + std::vector b(bytes); + for (auto & v : b) v = (uint8_t)(rnd() & 0xFF); + const size_t meta_off = block_bytes - 2; + for (size_t blk = 0; blk < bytes / block_bytes; ++blk) { + for (int h = 0; h < 2; ++h) { + uint8_t & m = b[blk * block_bytes + meta_off + h]; + // Scale index: a NARROW band, not the full UE4M3 range. Uniform + // random exponents span ~2^30 of dynamic range within one tensor, + // which no trained weight matrix does; against data like that any + // two accumulation orders disagree wildly and the comparison + // measures the fixture rather than the kernel. Real per-block + // scales in these artifacts sit within a couple of octaves of each + // other, so the band is chosen to match. The low bit stays random: + // it is the codebook-select flag and both settings must be + // exercised. + const uint8_t sel = (uint8_t)(m & 0x01); + m = (uint8_t)(0x3C + (rnd() % 9)) ; // ~2^-1 .. 2^1 + m = (uint8_t)((m & 0xFE) | sel); + } + } + return b; +} + +struct Err { double max_abs; double rms; }; + +Err compare(const std::vector & a, const std::vector & ref) { + double max_abs = 0.0, sq = 0.0; + for (size_t i = 0; i < ref.size(); ++i) { + const double d = std::fabs((double)a[i] - (double)ref[i]); + if (d > max_abs) max_abs = d; + sq += d * d; + } + return { max_abs, std::sqrt(sq / (double)ref.size()) }; +} + +// Run dst = src0^T * src1 through ggml on the CUDA/HIP backend, with the mix +// MMQ path forced on or off. Returns the dst rows. +bool run_ggml_mul_mat(ggml_backend_t backend, const MixType & mt, + const std::vector & blocks, + const std::vector & xh, + int in, int out, int n_experts, int ne11, + bool mmq, std::vector & out_rows) { + ggml_init_params ip{}; + ip.mem_size = ggml_tensor_overhead() * 8 + ggml_graph_overhead(); + ip.no_alloc = true; + ggml_context * ctx = ggml_init(ip); + if (!ctx) return false; + + ggml_tensor * w = ggml_new_tensor_3d(ctx, mt.type, in, out, n_experts); + ggml_tensor * x = ggml_new_tensor_3d(ctx, GGML_TYPE_F32, in, ne11, n_experts); + ggml_set_input(x); + ggml_tensor * y = ggml_mul_mat(ctx, w, x); + ggml_set_output(y); + ggml_cgraph * gf = ggml_new_graph(ctx); + ggml_build_forward_expand(gf, y); + + ggml_backend_buffer_t buf = ggml_backend_alloc_ctx_tensors(ctx, backend); + if (!buf) { ggml_free(ctx); return false; } + + ggml_backend_tensor_set(w, blocks.data(), 0, blocks.size()); + ggml_backend_tensor_set(x, xh.data(), 0, sizeof(float) * xh.size()); + + // Codebooks are keyed on the tensor's device pointer, so registration has + // to happen after allocation and be undone before the buffer is freed. + const int K = 8; + std::vector books((size_t)n_experts * 2 * K); + for (size_t i = 0; i < books.size(); ++i) { + const float v = 0.25f * (float)((int)(i % 7) - 3) + 0.05f * (float)(i / 16); + uint32_t bits; std::memcpy(&bits, &v, 4); + books[i] = (uint16_t)(bits >> 16); + } + std::vector modes(n_experts), rots(n_experts, 0); + for (int e = 0; e < n_experts; ++e) modes[(size_t)e] = (uint8_t)(e & 1); + + const size_t slice_bytes = blocks.size() / (size_t)n_experts; + mt.register_host(w->data, slice_bytes, n_experts, out, in, + books.data(), modes.data(), rots.data()); + + ggml_cuda_set_mix_mmq_enabled(mmq); + ggml_backend_graph_compute(backend, gf); + ggml_cuda_clear_mix_mmq_override(); + + out_rows.resize((size_t)out * ne11 * n_experts); + ggml_backend_tensor_get(y, out_rows.data(), 0, sizeof(float) * out_rows.size()); + + mt.unregister(w->data); + ggml_backend_buffer_free(buf); + ggml_free(ctx); + return true; +} + +// The MoE dispatch: ggml_mul_mat_id, which is what DeepSeek-V4's routed +// experts use. This is a DIFFERENT ladder from ggml_mul_mat above — a fused +// mix kernel for narrow batches, then mmvq, then MMQ, then a host-synchronised +// sort-and-dequantize fallback — and ggml_cuda_should_use_mmq gates the MMQ +// rung of it. Set DFLASH_MMID_TELEMETRY=1 to print the chosen path per call. +bool run_ggml_mul_mat_id(ggml_backend_t backend, const MixType & mt, + const std::vector & blocks, + const std::vector & xh, + const std::vector & ids_h, + int in, int out, int n_experts, int n_tokens, + int top_k, bool mmq, std::vector & out_rows) { + ggml_init_params ip{}; + ip.mem_size = ggml_tensor_overhead() * 8 + ggml_graph_overhead(); + ip.no_alloc = true; + ggml_context * ctx = ggml_init(ip); + if (!ctx) return false; + + // src0 [in, out, n_experts]; src1 [in, top_k, n_tokens]; ids [top_k, n_tokens] + ggml_tensor * w = ggml_new_tensor_3d(ctx, mt.type, in, out, n_experts); + ggml_tensor * x = ggml_new_tensor_3d(ctx, GGML_TYPE_F32, in, top_k, n_tokens); + ggml_tensor * ids = ggml_new_tensor_2d(ctx, GGML_TYPE_I32, top_k, n_tokens); + ggml_set_input(x); + ggml_set_input(ids); + ggml_tensor * y = ggml_mul_mat_id(ctx, w, x, ids); + ggml_set_output(y); + ggml_cgraph * gf = ggml_new_graph(ctx); + ggml_build_forward_expand(gf, y); + + ggml_backend_buffer_t buf = ggml_backend_alloc_ctx_tensors(ctx, backend); + if (!buf) { ggml_free(ctx); return false; } + + ggml_backend_tensor_set(w, blocks.data(), 0, blocks.size()); + ggml_backend_tensor_set(x, xh.data(), 0, sizeof(float) * xh.size()); + ggml_backend_tensor_set(ids, ids_h.data(), 0, sizeof(int32_t) * ids_h.size()); + + const int K = 8; + std::vector books((size_t)n_experts * 2 * K); + for (size_t i = 0; i < books.size(); ++i) { + const float v = 0.25f * (float)((int)(i % 7) - 3) + 0.05f * (float)(i / 16); + uint32_t bits; std::memcpy(&bits, &v, 4); + books[i] = (uint16_t)(bits >> 16); + } + std::vector modes(n_experts), rots(n_experts, 0); + for (int e = 0; e < n_experts; ++e) modes[(size_t)e] = (uint8_t)(e & 1); + + const size_t slice_bytes = blocks.size() / (size_t)n_experts; + mt.register_host(w->data, slice_bytes, n_experts, out, in, + books.data(), modes.data(), rots.data()); + + ggml_cuda_set_mix_mmq_enabled(mmq); + ggml_backend_graph_compute(backend, gf); + ggml_cuda_clear_mix_mmq_override(); + + out_rows.resize(ggml_nelements(y)); + ggml_backend_tensor_get(y, out_rows.data(), 0, sizeof(float) * out_rows.size()); + + mt.unregister(w->data); + ggml_backend_buffer_free(buf); + ggml_free(ctx); + return true; +} + +} // namespace + +int main() { + ggml_backend_t backend = ggml_backend_cuda_init(0); + if (!backend) { + std::fprintf(stderr, "SKIP: no GPU\n"); + return 77; + } + + const MixType types[] = { + { GGML_TYPE_Q3_1_ROCMFP3_MIX, "q3_1_rocmfp3_mix (105)", + ggml_cuda_rocmfp3_mix_register_host, ggml_cuda_rocmfp3_mix_unregister, + ggml_cuda_rocmfp3_mix_mul_mat_vec }, + { GGML_TYPE_Q2_1_ROCMFP2_MIX, "q2_1_rocmfp2_mix (106)", + ggml_cuda_rocmfp2_mix_register_host, ggml_cuda_rocmfp2_mix_unregister, + ggml_cuda_rocmfp2_mix_mul_mat_vec }, + }; + const int widths[] = { 1, 4, 16, 64 }; + // Dense only. A 3-D src0 (`n_experts > 1`) does NOT take the MMQ path in + // ggml — measured: the on and off runs came back bit-identical — so a + // multi-expert case here would compare a path against itself and pass + // while testing nothing. DS4's MoE weights go through `ggml_mul_mat_id`, + // a different dispatch, and are covered end-to-end against the real + // artifact instead. + const int expert_set[] = { 1 }; + + // The bar is NOT "as accurate as the dequant path". MMQ quantizes the + // activations to int8 where dequant+GEMM keeps them in f32, so MMQ is + // inherently coarser and measures ~6x the error of dequant on these types + // — stable across scale ranges, i.e. a property of the algorithms rather + // than of this fixture. Holding it to the dequant path's error would be + // demanding something MMQ cannot deliver, and llama.cpp ships MMQ by + // default for the K-quants on exactly this trade. + // + // What IS demanded is that the result is close to the reference in + // absolute terms: relative rms error below 1% of the output magnitude. + // The failure this exists to catch — a wrong codebook, a mis-strided + // tile, an ignored mode byte — puts the answer O(1) off, which is two + // orders of magnitude clear of this line. The MMQ/dequant ratio is + // reported for characterisation, not asserted on. + const double kRelRms = 0.01; + + for (const MixType & mt : types) { + const size_t block_bytes = ggml_type_size(mt.type); + const int qk = ggml_blck_size(mt.type); + for (int n_experts : expert_set) { + const int in = 256, out = 32; + const size_t nb = (size_t)(in / qk); + const size_t slice_bytes = (size_t)out * nb * block_bytes; + std::vector blocks = + make_blocks(slice_bytes * (size_t)n_experts, block_bytes); + + for (int ne11 : widths) { + std::vector xh((size_t)in * ne11 * n_experts); + for (auto & v : xh) v = 0.5f - (float)(rnd() % 1000) / 1000.0f; + + std::vector y_mmq, y_deq; + if (!run_ggml_mul_mat(backend, mt, blocks, xh, in, out, + n_experts, ne11, true, y_mmq) || + !run_ggml_mul_mat(backend, mt, blocks, xh, in, out, + n_experts, ne11, false, y_deq)) { + fail(std::string(mt.label) + ": graph run failed"); + continue; + } + + // Reference: the validated matvec kernel, one call per expert. + std::vector y_ref((size_t)out * ne11 * n_experts, 0.0f); + { + uint8_t * d_w = nullptr; float * d_x = nullptr; float * d_y = nullptr; + cudaMalloc(&d_w, blocks.size()); + cudaMemcpy(d_w, blocks.data(), blocks.size(), cudaMemcpyHostToDevice); + cudaMalloc(&d_x, sizeof(float) * xh.size()); + cudaMemcpy(d_x, xh.data(), sizeof(float) * xh.size(), cudaMemcpyHostToDevice); + cudaMalloc(&d_y, sizeof(float) * y_ref.size()); + cudaMemset(d_y, 0, sizeof(float) * y_ref.size()); + + const int K = 8; + std::vector books((size_t)n_experts * 2 * K); + for (size_t i = 0; i < books.size(); ++i) { + const float v = 0.25f * (float)((int)(i % 7) - 3) + 0.05f * (float)(i / 16); + uint32_t bits; std::memcpy(&bits, &v, 4); + books[i] = (uint16_t)(bits >> 16); + } + std::vector modes(n_experts), rots(n_experts, 0); + for (int e = 0; e < n_experts; ++e) modes[(size_t)e] = (uint8_t)(e & 1); + mt.register_host(d_w, slice_bytes, n_experts, out, in, + books.data(), modes.data(), rots.data()); + bool ok = true; + for (int e = 0; e < n_experts; ++e) { + ok &= mt.matvec(d_w + (size_t)e * slice_bytes, + d_x + (size_t)e * in * ne11, + d_y + (size_t)e * out * ne11, + in, out, ne11, in, out, nullptr); + } + cudaDeviceSynchronize(); + if (!ok) fail(std::string(mt.label) + ": reference matvec refused"); + cudaMemcpy(y_ref.data(), d_y, sizeof(float) * y_ref.size(), + cudaMemcpyDeviceToHost); + mt.unregister(d_w); + cudaFree(d_w); cudaFree(d_x); cudaFree(d_y); + } + + const Err e_mmq = compare(y_mmq, y_ref); + const Err e_deq = compare(y_deq, y_ref); + double mag = 0.0; + for (float v : y_ref) mag = std::max(mag, (double)std::fabs(v)); + std::printf("%-24s experts=%d ne11=%-3d |ref|max %.4g | MMQ rms %.4g " + "(%.3f%% of |ref|) | dequant rms %.4g | ratio %.1fx\n", + mt.label, n_experts, ne11, mag, + e_mmq.rms, mag > 0 ? 100.0 * e_mmq.rms / mag : 0.0, + e_deq.rms, e_deq.rms > 0 ? e_mmq.rms / e_deq.rms : 0.0); + + // Guard against a vacuous pass: if the reference is all zeros + // (a decode that produced nothing), every error is 0 and the + // comparison means nothing. + double ref_mag = 0.0; + for (float v : y_ref) ref_mag = std::max(ref_mag, (double)std::fabs(v)); + if (ref_mag < 1e-6) { + fail(std::string(mt.label) + ": reference output is all-zero, " + "the comparison would pass vacuously"); + continue; + } + // Did MMQ actually engage? For ne11 > 1 the two runs take + // different kernels and cannot agree bit-for-bit; if they do, + // the toggle did nothing and this row is comparing a path + // against itself. That is exactly how the first version of + // this test "passed" every multi-expert case. + if (ne11 > 1 && y_mmq == y_deq) { + fail(std::string(mt.label) + " experts=" + + std::to_string(n_experts) + " ne11=" + + std::to_string(ne11) + ": MMQ and dequant runs are " + "bit-identical, so MMQ never engaged and this case " + "tests nothing"); + continue; + } + if (e_mmq.rms > kRelRms * ref_mag) { + fail(std::string(mt.label) + " experts=" + std::to_string(n_experts) + + " ne11=" + std::to_string(ne11) + ": MMQ rms " + + std::to_string(e_mmq.rms) + " is " + + std::to_string(100.0 * e_mmq.rms / ref_mag) + + "% of |ref| max " + std::to_string(ref_mag) + + ", over the 1% bar"); + } + } + } + } + + // ---- MoE: ggml_mul_mat_id, the path DeepSeek-V4's routed experts take ---- + // + // The dense loop above deliberately uses n_experts == 1: a 3-D src0 through + // ggml_mul_mat does not reach MMQ, so a multi-expert case THERE would + // compare a path against itself. The real MoE dispatch is ggml_mul_mat_id, + // which has its own ladder (fused mix kernel for narrow batches -> mmvq -> + // MMQ -> host-synchronised sort fallback) and which ggml_cuda_should_use_mmq + // gates. Leaving it untested is what let "expert counts > 1 are covered" be + // claimed when they were not. + // + // The cross-check here is MMQ against dequantize+GEMM through the same + // mul_mat_id call: two independent implementations of one operation. The + // widths straddle the ladder deliberately, and the report says for each + // width whether the toggle changed the result at all — which is the only + // way to see, from outside, where MMQ actually engages for MoE. + { + const int n_experts = 8, top_k = 2; + const int in = 256, out = 32; + const int token_widths[] = { 1, 4, 64, 512 }; + for (const MixType & mt : types) { + const size_t block_bytes = ggml_type_size(mt.type); + const int qk = ggml_blck_size(mt.type); + const size_t nb = (size_t)(in / qk); + const size_t slice_bytes = (size_t)out * nb * block_bytes; + std::vector blocks = + make_blocks(slice_bytes * (size_t)n_experts, block_bytes); + + for (int n_tokens : token_widths) { + std::vector xh((size_t)in * top_k * n_tokens); + for (auto & v : xh) v = 0.5f - (float)(rnd() % 1000) / 1000.0f; + // A token's top-k picks must be DISTINCT: the runtime's sort + // fallback records one entry per (expert, token) and breaks on + // the first match, so a duplicated expert makes it emit fewer + // rows than ne12*n_expert_used and trip its own assert. Real + // top-k routing cannot select an expert twice; sampling with + // replacement here produced an abort that looked like a kernel + // bug and was bad fixture data. + std::vector ids((size_t)top_k * n_tokens); + { + std::vector pool(n_experts); + for (int e = 0; e < n_experts; ++e) pool[(size_t)e] = e; + for (int t = 0; t < n_tokens; ++t) { + for (int k = 0; k < top_k; ++k) { + const int j = k + (int)(rnd() % (unsigned)(n_experts - k)); + std::swap(pool[(size_t)k], pool[(size_t)j]); + ids[(size_t)t * top_k + k] = pool[(size_t)k]; + } + } + } + + std::vector y_mmq, y_deq; + if (!run_ggml_mul_mat_id(backend, mt, blocks, xh, ids, in, out, + n_experts, n_tokens, top_k, true, y_mmq) || + !run_ggml_mul_mat_id(backend, mt, blocks, xh, ids, in, out, + n_experts, n_tokens, top_k, false, y_deq)) { + fail(std::string(mt.label) + " mul_mat_id: graph run failed"); + continue; + } + + // Timed repeat: which is actually faster at this width, the + // MMQ path or whatever declining it lands on? For mul_mat_id + // that fallback is the host-synchronised sort, not a dense GEMM. + double t_mmq = 0.0, t_deq = 0.0; + { + std::vector scratch; + const int reps = 5; + auto clk = []() { + return std::chrono::duration( + std::chrono::steady_clock::now().time_since_epoch()).count(); + }; + run_ggml_mul_mat_id(backend, mt, blocks, xh, ids, in, out, + n_experts, n_tokens, top_k, true, scratch); + const double a0 = clk(); + for (int r = 0; r < reps; ++r) + run_ggml_mul_mat_id(backend, mt, blocks, xh, ids, in, out, + n_experts, n_tokens, top_k, true, scratch); + t_mmq = (clk() - a0) / reps; + run_ggml_mul_mat_id(backend, mt, blocks, xh, ids, in, out, + n_experts, n_tokens, top_k, false, scratch); + const double b0 = clk(); + for (int r = 0; r < reps; ++r) + run_ggml_mul_mat_id(backend, mt, blocks, xh, ids, in, out, + n_experts, n_tokens, top_k, false, scratch); + t_deq = (clk() - b0) / reps; + } + + const Err e = compare(y_mmq, y_deq); + double mag = 0.0; + for (float v : y_deq) mag = std::max(mag, (double)std::fabs(v)); + const bool engaged = (y_mmq != y_deq); + std::printf("%-24s mul_mat_id experts=%d top_k=%d tokens=%-4d | " + "MMQ %s | rms vs dequant %.5g (%.3f%% of |ref| %.4g)\n", + mt.label, n_experts, top_k, n_tokens, + engaged ? "engaged " : "NOT engaged", e.rms, + mag > 0 ? 100.0 * e.rms / mag : 0.0, mag); + std::printf("%-24s timing: mmq-on %.3f ms | mmq-off %.3f ms | %s\n", + mt.label, t_mmq * 1e3, t_deq * 1e3, + t_mmq < t_deq ? "MMQ-ON faster" : "MMQ-OFF faster"); + + // Only assert accuracy where the toggle actually changed the + // kernel. Where it did not, there is nothing to compare and + // saying "pass" would be the vacuous result this guards against. + if (engaged && e.rms > kRelRms * mag) { + fail(std::string(mt.label) + " mul_mat_id tokens=" + + std::to_string(n_tokens) + ": MMQ rms " + + std::to_string(e.rms) + " is " + + std::to_string(mag > 0 ? 100.0 * e.rms / mag : 0.0) + + "% of |ref|, over the 1% bar"); + } + } + } + } + + ggml_backend_free(backend); + if (g_fails == 0) { std::printf("test_rocmfp_mix_mmq: OK\n"); return 0; } + std::fprintf(stderr, "test_rocmfp_mix_mmq: %d failure(s)\n", g_fails); + return 1; +}