cuda: fuse MoE weighted expert reduction (#25952)

* cuda : fuse MoE weighted reduction (mul + view + add)

The MoE combine tail currently writes weighted expert outputs to
global memory before reducing them. That intermediate global-memory
traffic is the main cost. The production baseline generally runs two
physical fused kernels; this path runs one.

This change matches the full expert-weighting plus ordered-reduction
subgraph and replaces it with one weighted-reduction kernel.

Supported graphs:
- unscaled: experts * router_weights
- scaled:   (experts * expert_scale) * router_weights

k = 2..15 is handled by one runtime-k kernel.

Matching is structural: op sequence, shapes, strides, expert views,
and the left-to-right ADD chain. The fused kernel keeps that same
reduction order. Results are not claimed bit-identical; CUDA FP32
contraction can change rounding slightly.

Allocator integration uses add_alloc_dep from the graph-optimizer
API so experts, router weights, and optional expert scales stay live
until the fused destination is written. Memory ranges are rechecked
before the fused kernel runs.

Unrecognized or unsafe graphs are left alone and keep the existing
per-op path. Set GGML_CUDA_MOE_WEIGHTED_REDUCTION=0 to disable the
fusion.

test-backend-ops covers scaled/unscaled, aligned/unaligned, and
representative values across k=2..15, plus a k=16 case that must
stay on the per-op path.

* Pruned the test matrix from 15 to 6

* Addressed the aman and olivers review comments

ggml/src/ggml-cuda/moe-weighted-reduction.cu

@@ -0,0 +1,65 @@
+#include "moe-weighted-reduction.cuh"
+
+static __global__ void moe_weighted_reduction_f32(const float * __restrict__ experts,
+                                                  const float * __restrict__ expert_scale,
+                                                  const float * __restrict__ weights,
+                                                  float * __restrict__ dst,
+                                                  const int64_t n_embd,
+                                                  const int     n_expert_used) {
+    const int64_t token = blockIdx.x;
+    const int64_t col   = (int64_t) blockIdx.y * blockDim.x + threadIdx.x;
+    if (col >= n_embd) {
+        return;
+    }
+
+    const uint64_t first_row   = (uint64_t) token * n_expert_used;
+    const float    first_scale = expert_scale != nullptr ? expert_scale[first_row] : 1.0f;
+    float          sum         = (experts[first_row * n_embd + col] * first_scale) * weights[first_row];
+
+    for (int expert = 1; expert < n_expert_used; ++expert) {
+        const uint64_t row   = first_row + expert;
+        const float   scale = expert_scale != nullptr ? expert_scale[row] : 1.0f;
+        sum += (experts[row * n_embd + col] * scale) * weights[row];
+    }
+    dst[token * n_embd + col] = sum;
+}
+
+static void launch_moe_weighted_reduction(const float * experts,
+                                          const float * expert_scale,
+                                          const float * weights,
+                                          float *       dst,
+                                          int64_t       n_embd,
+                                          int64_t       n_tokens,
+                                          int           n_expert_used,
+                                          cudaStream_t  stream) {
+    constexpr int threads = 256;
+    const dim3 blocks(n_tokens, (n_embd + threads - 1) / threads, 1);
+    moe_weighted_reduction_f32
+        <<<blocks, threads, 0, stream>>>(experts, expert_scale, weights, dst, n_embd, n_expert_used);
+}
+
+void ggml_cuda_op_moe_weighted_reduction(ggml_backend_cuda_context & ctx,
+                                         const ggml_tensor *         experts,
+                                         const ggml_tensor *         expert_scale,
+                                         const ggml_tensor *         weights,
+                                         ggml_tensor *               dst) {
+    GGML_ASSERT(experts->type == GGML_TYPE_F32);
+    GGML_ASSERT(weights->type == GGML_TYPE_F32);
+    GGML_ASSERT(expert_scale == nullptr || expert_scale->type == GGML_TYPE_F32);
+    GGML_ASSERT(dst->type == GGML_TYPE_F32);
+    GGML_ASSERT(ggml_is_contiguous(experts));
+    GGML_ASSERT(ggml_is_contiguous(weights));
+    GGML_ASSERT(expert_scale == nullptr || ggml_is_contiguous(expert_scale));
+    GGML_ASSERT(ggml_is_contiguous(dst));
+
+    const int64_t n_embd        = experts->ne[0];
+    const int64_t n_expert_used = experts->ne[1];
+    const int64_t n_tokens      = experts->ne[2] * experts->ne[3];
+    cudaStream_t  stream        = ctx.stream();
+
+    launch_moe_weighted_reduction((const float *) experts->data,
+                                  expert_scale ? (const float *) expert_scale->data : nullptr,
+                                  (const float *) weights->data,
+                                  (float *) dst->data, n_embd, n_tokens, (int) n_expert_used, stream);
+    CUDA_CHECK(cudaGetLastError());
+}