ggml : add fused SINKHORN_NORM op (Sinkhorn-Knopp normalization) + use in DeepSeek-V4

## Summary

This adds a new ggml op, `GGML_OP_SINKHORN_NORM` (`ggml_sinkhorn_norm`), that performs a
fused **Sinkhorn-Knopp doubly-stochastic normalization** on small square matrices, and wires
it into the DeepSeek-V4 (mHC / manifold-constrained hyper-connections) graph.

DeepSeek-V4 makes its per-token, per-layer hyper-connection mixing matrix (`[hc, hc]`,
with `hc = 4`) doubly-stochastic by running Sinkhorn-Knopp (softmax, then alternating
row/column normalization) for `sinkhorn_iterations` (20 in the released model). Today this
is expressed as a decomposition of primitive ggml ops (`soft_max` + repeated
`sum_rows` / `add` / `div` / `cont(permute)` transpose), which expands to **~137 ops per
layer**. On the tiny 4x4 slices this is almost entirely kernel-launch and memory-round-trip
overhead with negligible arithmetic.

The new op computes the entire normalization in a single kernel (per layer), keeping the
matrix in registers across all iterations.

## Motivation / profiling

Profiling DeepSeek-V4-Flash (UD-Q4_K_XL, 284B, 43 layers) inference on an AMD MI250X
(gfx90a, ROCm/HIP) with `rocprofv3 --kernel-trace` showed the decomposed Sinkhorn dominating
GPU time:

- `k_bin_bcast` (broadcast div/add): ~31% of GPU kernel time, ~217k dispatches
- `reduce_rows_f32` (row sums): ~8%, ~67k dispatches
- `cpy_scalar_transpose`: ~4.5%, ~34k dispatches

Predicted vs. observed dispatch counts confirmed these are almost entirely the Sinkhorn
decomposition (e.g. transpose count matched to within <1%). In aggregate the normalization
was **~30% of total GPU kernel time** across ~250k tiny dispatches.

## Results (DeepSeek-V4-Flash, 4x MI250X GCDs, llama-bench, back-to-back A/B)

| variant                          | prefill (pp64) | decode (tg16) |
|----------------------------------|---------------:|--------------:|
| decomposed (baseline)            |   111.5 t/s    |   10.3 t/s    |
| fused SINKHORN_NORM (this PR)     |   119.8 t/s    |   16.5 t/s    |

**Decode +60.9%, prefill +7.4%** end-to-end on the full model.

## Implementation

- New op `GGML_OP_SINKHORN_NORM` + constructor
  `ggml_sinkhorn_norm(ctx, a, n_iters, eps)` (op params: `n_iters` int32, `eps` f32).
  Operates independently on each `[n, n]` slice (`ne0 == ne1 == n`): softmax over `ne0`,
  add `eps`, one column normalization, then `n_iters - 1` alternating {row, column}
  normalizations (each `+ eps`). The interface is general (arbitrary square `n`, configurable
  `n_iters` / `eps`); DeepSeek-V4 is simply the first consumer.
- **CPU** reference implementation (`ggml-cpu/ops.cpp`), parallelized over slices.
- **CUDA/HIP** kernel (`ggml-cuda/sinkhorn-norm.cu`):
  - warp-cooperative kernel for power-of-two `n` with `n*n <= warp_size`: one lane per matrix
    element, fully coalesced global load/store, segmented `__shfl_xor` butterfly reductions
    for the row/column sums;
  - a generic one-thread-per-slice fallback otherwise.
- **Other backends** (Metal/Vulkan/SYCL/CANN/OpenCL): `supports_op` returns false for this
  op (default), so the graph scheduler falls back to the CPU implementation. No crashes; a
  dedicated kernel can be added later if desired.
- **Tests**: `test-backend-ops` cases for `SINKHORN_NORM` (n=4, tokens {1,7,128,1024},
  iters {1,2,20}) plus perf cases; CPU-vs-CUDA all pass.
- DeepSeek-V4 (`src/models/deepseek4.cpp`) `build_hc_sinkhorn` now calls the fused op.

## Precedent

ggml already contains architecture-specific fused ops in the core op enum, e.g.
`GGML_OP_SSM_CONV` / `GGML_OP_SSM_SCAN` (Mamba), `GGML_OP_RWKV_WKV6` / `GGML_OP_RWKV_WKV7` /
`GGML_OP_GATED_LINEAR_ATTN` (RWKV). `SINKHORN_NORM` follows the same pattern for DeepSeek-V4,
but is defined as a general normalization primitive rather than a model-specific op.

## Test plan

- [ ] `test-backend-ops test -o SINKHORN_NORM` passes on CPU and CUDA/HIP (12/12 here).
- [ ] `test-backend-ops perf -o SINKHORN_NORM` runs.
- [ ] Full `test-backend-ops` shows no regressions (new op cleanly inserted into the enum).
- [ ] DeepSeek-V4 output unchanged vs. the decomposed graph (validated via the CPU reference).
- [ ] Confirm graceful CPU fallback on a non-CUDA backend (e.g. Metal/Vulkan).

## Notes for reviewers

- Open question / happy to adjust: keep this as a public `GGML_OP_SINKHORN_NORM`, or fold it
  into an internal fused helper for DeepSeek-V4 only? Went with a general op given the
  Mamba/RWKV precedent, but will follow maintainer preference.
- Warp kernel assumes `n` is a power of two and `n*n <= warp_size` (true for `hc = 4`);
  the generic kernel covers everything else.

Author: jadenmach2

ggml/include/ggml.h

@@ -509,6 +509,7 @@ extern "C" {
         GGML_OP_RMS_NORM_BACK,
         GGML_OP_GROUP_NORM,
         GGML_OP_L2_NORM,
+        GGML_OP_SINKHORN_NORM,
 
         GGML_OP_MUL_MAT,
         GGML_OP_MUL_MAT_ID,
@@ -1406,6 +1407,17 @@ extern "C" {
             struct ggml_tensor  * a,
             float                 eps);
 
+    // fused Sinkhorn-Knopp normalization for DeepSeek-V4 mHC hyper-connections.
+    // operates independently on each [n, n] slice (ne0 == ne1 == n) of the input:
+    //   softmax over ne0, add eps, then a row normalization followed by
+    //   (n_iters - 1) alternating {column, row} normalizations (each + eps).
+    // replaces the ~137-op decomposed graph with a single kernel per call.
+    GGML_API struct ggml_tensor * ggml_sinkhorn_norm(
+            struct ggml_context * ctx,
+            struct ggml_tensor  * a,
+            int                   n_iters,
+            float                 eps);
+
     // a - x
     // b - dy
     GGML_API struct ggml_tensor * ggml_rms_norm_back(

ggml/src/ggml-cpu/ggml-cpu.c

@@ -1833,6 +1833,10 @@ static void ggml_compute_forward(struct ggml_compute_params * params, struct ggm
             {
                 ggml_compute_forward_l2_norm(params, tensor);
             } break;
+        case GGML_OP_SINKHORN_NORM:
+            {
+                ggml_compute_forward_sinkhorn_norm(params, tensor);
+            } break;
         case GGML_OP_MUL_MAT:
             {
                 ggml_compute_forward_mul_mat(params, tensor);
@@ -2306,6 +2310,7 @@ static int ggml_get_n_tasks(struct ggml_tensor * node, int n_threads) {
         case GGML_OP_RMS_NORM:
         case GGML_OP_RMS_NORM_BACK:
         case GGML_OP_L2_NORM:
+        case GGML_OP_SINKHORN_NORM:
         case GGML_OP_GROUP_NORM:
         case GGML_OP_CONCAT:
         case GGML_OP_MUL_MAT:

ggml/src/ggml-cpu/ops.cpp

@@ -4237,6 +4237,131 @@ void ggml_compute_forward_l2_norm(
     }
 }
 
+// ggml_compute_forward_sinkhorn_norm
+
+#define GGML_SINKHORN_NORM_MAX_N 32
+
+static void ggml_compute_forward_sinkhorn_norm_f32(
+    const ggml_compute_params * params,
+    ggml_tensor * dst) {
+
+    const ggml_tensor * src0 = dst->src[0];
+
+    GGML_ASSERT(ggml_are_same_shape(src0, dst));
+    GGML_ASSERT(ggml_is_contiguous(src0));
+
+    const int ith = params->ith;
+    const int nth = params->nth;
+
+    GGML_TENSOR_UNARY_OP_LOCALS
+
+    const int32_t n_iters = ggml_get_op_params_i32(dst, 0);
+    const float   eps     = ggml_get_op_params_f32(dst, 1);
+
+    const int64_t n = ne00; // square [n, n] slices, ne00 == ne01
+    GGML_ASSERT(ne01 == n);
+    GGML_ASSERT(n <= GGML_SINKHORN_NORM_MAX_N);
+    GGML_ASSERT(n_iters >= 1);
+
+    const int64_t n_slices = ne02 * ne03;
+
+    float m[GGML_SINKHORN_NORM_MAX_N * GGML_SINKHORN_NORM_MAX_N];
+
+    // one square slice (spanning ne0 x ne1) per token; parallelize across slices
+    for (int64_t s = ith; s < n_slices; s += nth) {
+        const int64_t i02 = s % ne02;
+        const int64_t i03 = s / ne02;
+
+        const float * x = (const float *) ((const char *) src0->data + i02*nb02 + i03*nb03);
+        float *       y = (float *)       ((      char *) dst->data  + i02*nb2  + i03*nb3);
+
+        // load slice: m[b*n + a], a is ne0 (fastest, contiguous)
+        for (int64_t b = 0; b < n; ++b) {
+            for (int64_t a = 0; a < n; ++a) {
+                m[b*n + a] = x[b*n + a];
+            }
+        }
+
+        // softmax over ne0 (a) for each b (matches ggml_soft_max, with max subtraction)
+        for (int64_t b = 0; b < n; ++b) {
+            float mx = -INFINITY;
+            for (int64_t a = 0; a < n; ++a) {
+                mx = fmaxf(mx, m[b*n + a]);
+            }
+            float sum = 0.0f;
+            for (int64_t a = 0; a < n; ++a) {
+                const float e = expf(m[b*n + a] - mx);
+                m[b*n + a] = e;
+                sum += e;
+            }
+            for (int64_t a = 0; a < n; ++a) {
+                m[b*n + a] /= sum;
+            }
+        }
+
+        // add eps
+        for (int64_t i = 0; i < n*n; ++i) {
+            m[i] += eps;
+        }
+
+        // norm_cols: divide each element by (sum over b of its row a) + eps
+        auto norm_cols = [&]() {
+            for (int64_t a = 0; a < n; ++a) {
+                float r = 0.0f;
+                for (int64_t b = 0; b < n; ++b) {
+                    r += m[b*n + a];
+                }
+                r += eps;
+                for (int64_t b = 0; b < n; ++b) {
+                    m[b*n + a] /= r;
+                }
+            }
+        };
+
+        // norm_rows: divide each element by (sum over a of its column b) + eps
+        auto norm_rows = [&]() {
+            for (int64_t b = 0; b < n; ++b) {
+                float c = 0.0f;
+                for (int64_t a = 0; a < n; ++a) {
+                    c += m[b*n + a];
+                }
+                c += eps;
+                for (int64_t a = 0; a < n; ++a) {
+                    m[b*n + a] /= c;
+                }
+            }
+        };
+
+        norm_cols();
+        for (int32_t it = 1; it < n_iters; ++it) {
+            norm_rows();
+            norm_cols();
+        }
+
+        for (int64_t i = 0; i < n*n; ++i) {
+            y[i] = m[i];
+        }
+    }
+}
+
+void ggml_compute_forward_sinkhorn_norm(
+    const ggml_compute_params * params,
+    ggml_tensor * dst) {
+
+    const ggml_tensor * src0 = dst->src[0];
+
+    switch (src0->type) {
+        case GGML_TYPE_F32:
+            {
+                ggml_compute_forward_sinkhorn_norm_f32(params, dst);
+            } break;
+        default:
+            {
+                GGML_ABORT("fatal error");
+            }
+    }
+}
+
 // ggml_compute_forward_out_prod
 
 static void ggml_compute_forward_out_prod_f32(

ggml/src/ggml-cpu/ops.h

@@ -48,6 +48,7 @@ void ggml_compute_forward_rms_norm_mul_fused(const struct ggml_compute_params *
 void ggml_compute_forward_rms_norm_back(const struct ggml_compute_params * params, struct ggml_tensor * dst);
 void ggml_compute_forward_group_norm(const struct ggml_compute_params * params, struct ggml_tensor * dst);
 void ggml_compute_forward_l2_norm(const struct ggml_compute_params * params, struct ggml_tensor * dst);
+void ggml_compute_forward_sinkhorn_norm(const struct ggml_compute_params * params, struct ggml_tensor * dst);
 void ggml_compute_forward_out_prod(const struct ggml_compute_params * params, struct ggml_tensor * dst);
 void ggml_compute_forward_scale(const struct ggml_compute_params * params, struct ggml_tensor * dst);
 void ggml_compute_forward_set(const struct ggml_compute_params * params, struct ggml_tensor * dst);

ggml/src/ggml-cuda/ggml-cuda.cu

@@ -33,6 +33,7 @@
 #include "ggml-cuda/mmvf.cuh"
 #include "ggml-cuda/mmvq.cuh"
 #include "ggml-cuda/norm.cuh"
+#include "ggml-cuda/sinkhorn-norm.cuh"
 #include "ggml-cuda/opt-step-adamw.cuh"
 #include "ggml-cuda/opt-step-sgd.cuh"
 #include "ggml-cuda/out-prod.cuh"
@@ -2083,6 +2084,9 @@ static bool ggml_cuda_compute_forward(ggml_backend_cuda_context & ctx, struct gg
         case GGML_OP_GROUP_NORM:
             ggml_cuda_op_group_norm(ctx, dst);
             break;
+        case GGML_OP_SINKHORN_NORM:
+            ggml_cuda_op_sinkhorn_norm(ctx, dst);
+            break;
         case GGML_OP_L2_NORM:
             ggml_cuda_op_l2_norm(ctx, dst);
             break;
@@ -4841,6 +4845,11 @@ static bool ggml_backend_cuda_device_supports_op(ggml_backend_dev_t dev, const g
         case GGML_OP_SILU_BACK:
             return ggml_is_contiguous(op->src[0]) && op->src[0]->type == GGML_TYPE_F32;
             break;
+        case GGML_OP_SINKHORN_NORM:
+            return op->src[0]->type == GGML_TYPE_F32 &&
+                   ggml_is_contiguous(op->src[0]) &&
+                   op->src[0]->ne[0] == op->src[0]->ne[1] &&
+                   op->src[0]->ne[0] <= 8;
         case GGML_OP_NORM:
         case GGML_OP_RMS_NORM:
         case GGML_OP_L2_NORM: