[CUDA] Add decode (M=1) GEMV fast path to MatMul

[Tianlei Wu (tianleiwu)]

Description

Adds a decode-oriented (M == 1) GEMV fast path to the CUDA MatMul op. At decode a plain
fp16/bf16 MatMul with a constant weight is a memory-bound GEMV; cuBLAS dispatches a
split-K path (a dot_kernel followed by a separate reduce_1Block_kernel) — two launches
for one tiny GEMV. A single custom kernel that prepacks the constant weight (transposed for
coalescing) and assigns one block per output column is ~2× faster.

Motivation

Profiling decode of a hybrid linear-attention + MoE model (Qwen3-Next / Qwen3.6) showed the
MoE router-gate MatMul nodes ([1, 2048] × [2048, 256] and a [2048, 1] projection,
fp16) spending ~350 µs/token across the model. They are plain MatMul ops that ORT routes
to cuBLAS, which selects a 496-way split-K (dot_kernel) + separate reduce_1Block_kernel
— two kernels for one small M = 1 GEMV.

A microbench of [1, 2048] × [2048, 256] fp16 showed:

Variant	µs/call	vs cuBLAS
cuBLAS Hgemm (split-K dot + reduce)	7.11	1.0×
custom GEMV, weight row-major [K, N] (strided)	7.45	0.96×
custom GEMV, weight transposed [N, K] (coalesced)	3.65	1.95×

The row-major variant ties cuBLAS — the coalesced transposed layout is what wins, so the
weight is transposed once during PrePack.

Key Changes

File	Change
core/providers/cuda/math/matmul_gemv.{h,cu}	New GemvM1Kernel (one block per output column, splits K across the block, single block-reduce, fp32 accumulation) + a tiled TransposeForGemv helper. Instantiated for fp16/bf16 only.
core/providers/cuda/math/matmul.{h,cc}	PrePack transposes an eligible constant B into an [N, K] buffer; ComputeInternal dispatches the GEMV at M == 1.
test/providers/cpu/math/matmul_test.cc	MatMul_Float16_GemvDecode parity test.

Eligibility (all required): input is the constant weight B (initializer, input 1),
dtype fp16/bf16, used as a plain 2-D [K, N] matrix (no transB/transBatchB),
N ≤ 1024, K ≥ 256. is_packed is left false so the original B tensor stays
available — at M > 1 (and for transposed A, batched GEMM, or non-eligible shapes) the
kernel falls through to the existing cuBLAS path unchanged.

Performance

End-to-end on Qwen3.6-35B-A3B (H200, INT4, single-sequence decode, CUDA graph on):

• The router-gate cuBLAS kernels (nvjet 207 + dot_kernel 73 + reduce_1Block 72 =
  ~352 µs/token) are replaced by a single GemvM1Kernel (~192 µs/token).
• Total GPU busy 5.25 → 5.10 ms/token.
• Decode throughput 154.8 → 158.3 tok/s (+2.3%).

Testing

MatMul_Float16_GemvDecode exercises both paths (constant-B fast path and non-constant-B
cuBLAS path) against a reference; fp32 accumulation matches cuBLAS's fp32-compute fp16 GEMM.

./onnxruntime_provider_test --gtest_filter='MathOpTest.MatMul*:*FusedMatMul*'

All 36 MatMul / FusedMatMul tests pass.

Motivation and Context

Decode-throughput optimization for fp16/bf16 models with small constant-weight projections
(e.g. MoE router gates). No public API change; numerics are preserved and all non-eligible
shapes keep the existing cuBLAS behavior.

[github-actions]

You can commit the suggested changes from lintrunner.