Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
326 changes: 93 additions & 233 deletions README.md

Large diffs are not rendered by default.

68 changes: 68 additions & 0 deletions scripts/bench-129-ab.ps1
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
# Issue #129 A/B: GPU-SLRU MoE prefill (SHARPI_CPU_MOE=0) — this branch vs master.
# Forces the on-GPU routed-expert path that #129's fused weighted-scatter-reduce kernel
# optimizes (the default auto-routes MoE to CPU, so no standard bench row exercises it).
# Confirms Backend: CUDA hybrid + GPU-SLRU in EACH run before trusting the number
# (a prior $args/PowerShell collision silently fell back to CPU).
param(
[string]$Model = "E:\models\Qwen3.6-35B-A3B-UD-Q4_K_M.gguf",
[int]$NTokens = 8
)

# Identical ~1K-token prompt to bench-allrows-1k.ps1.
$para = @(
"Modern large language model inference is dominated by memory bandwidth rather than raw compute.",
"Each decode step streams the full weight matrix for every layer, so quantization formats like Q4_K and Q5_K trade a small accuracy loss for a large reduction in bytes moved per token.",
"Mixture-of-experts models complicate this: only a handful of the hundreds of experts fire per token, and which fire varies token to token, defeating simple weight caching.",
"Gated DeltaNet layers replace quadratic attention with a linear recurrent state update, bounding per-token cost as context grows, at the price of a strictly sequential scan.",
"Hybrid placement keeps the attention and recurrent trunk on the accelerator while streaming routed-expert weights from host memory, overlapping the two so neither stalls."
) -join " "
$sb = [System.Text.StringBuilder]::new()
[void]$sb.Append("Read the following engineering notes and then write a concise technical summary.`n`n")
for ($i = 1; $i -le 6; $i++) { [void]$sb.Append("Section $i. $para`n`n") }
[void]$sb.Append("Summarize the main performance trade-offs across the sections above.")
$prompt = $sb.ToString()

$dlls = [ordered]@{
"master" = "C:\p\sharpi-master\src\SharpInference.Cli\bin\Release\net10.0\sharpi-cli.dll"
"branch129" = "C:\p\sharpi\src\SharpInference.Cli\bin\Release\net10.0\sharpi-cli.dll"
}

$env:SHARPI_CPU_MOE = "0" # force GPU-SLRU on-GPU routed experts (#129 path)
$dotnet = (Get-Command dotnet).Source
$outDir = "C:\p\sharpi\tools\bench"
if (-not (Test-Path $outDir)) { New-Item -ItemType Directory -Path $outDir | Out-Null }

function Run([string]$label, [string]$dll, [bool]$warm) {
$args = @($dll, "-m", $Model, "-p", $prompt, "--temp", "0", "-n", "$NTokens",
"-g", "-1", "--backend", "cuda", "--single-turn", "--verbose-prompt")
$psi = New-Object System.Diagnostics.ProcessStartInfo
$psi.FileName = $dotnet
$psi.Arguments = ($args | ForEach-Object { if ($_ -match '\s') { "`"$_`"" } else { $_ } }) -join ' '
$psi.RedirectStandardOutput = $true; $psi.RedirectStandardError = $true; $psi.UseShellExecute = $false
$p = [System.Diagnostics.Process]::Start($psi)
$so = $p.StandardOutput.ReadToEndAsync(); $se = $p.StandardError.ReadToEndAsync()
$p.WaitForExit(900*1000) | Out-Null
$so.Wait(3000) | Out-Null; $se.Wait(3000) | Out-Null
$out = $so.Result + "`n" + $se.Result
if ($warm) { return }
Set-Content -Path (Join-Path $outDir "129ab-$label.txt") -Value $out
# Backend confirmation (contract: must be CUDA hybrid + GPU-SLRU, NOT CPU fallback).
$slru = ($out | Select-String -Pattern "SLRU expert cache").Matches.Count
$cpuMoe = ($out | Select-String -Pattern "MoE.*on CPU|auto-routed to CPU|Dense FFN mode").Matches.Count
$backendLine = ($out -split "`n" | Where-Object { $_ -match "Backend|SLRU expert cache|GDN:|MoE:" } | Select-Object -First 4) -join " | "
$prefill = "n/a"
if ($out -match 'Prefill:\s+(\d+)\s+tokens,\s+([\d\.]+)\s+t/s') { $prefill = "$($matches[1]) tok @ $($matches[2]) t/s" }
$decode = "n/a"
if ($out -match 'Decode:\s+(\d+)\s+tokens,\s+([\d\.]+)\s+t/s') { $decode = "$($matches[2]) t/s" }
[PSCustomObject]@{ Label=$label; Prefill=$prefill; Decode=$decode; SLRUlines=$slru; CpuMoeHits=$cpuMoe; Backend=$backendLine }
}

Write-Host "=== Warming OS page cache for $Model ===" -ForegroundColor DarkGray
Run "warm" $dlls["branch129"] $true

$results = @()
foreach ($kv in $dlls.GetEnumerator()) {
Write-Host "=== Running $($kv.Key) (SHARPI_CPU_MOE=0, GPU-SLRU) ===" -ForegroundColor Cyan
$results += Run $kv.Key $kv.Value $false
}
$results | Format-List
71 changes: 71 additions & 0 deletions src/SharpInference.Cuda/CudaBackend.cs
Original file line number Diff line number Diff line change
Expand Up @@ -163,6 +163,12 @@ public sealed unsafe class CudaBackend : IComputeBackend, IImageOpsBackend, IDis
private nint _clearF32Kernel;
private nint _quantizeQ81Kernel;
private nint _bwBaselineKernel;
// Issue #129: batched GPU-SLRU MoE host-loop replacements. scale_rows applies a
// per-row scalar (shared-expert sigmoid gate) over [rows × cols] in one launch;
// moe_weighted_reduce does the whole Phase-3 top-k weighted scatter-reduce +
// shared add over all N tokens in one launch (replacing ~N·(na+2) tiny ops).
private nint _scaleRowsKernel;
private nint _moeWeightedReduceKernel;

// TurboQuant KV-cache compression kernels.
private nint _tqRotateQueryKernel;
Expand Down Expand Up @@ -3453,6 +3459,7 @@ private void ForceEagerJit()
_attentionKernel, _attentionBf16Kernel, _attentionSwaKernel,
_geluTanhMulKernel, _softcapKernel,
_clearF32Kernel, _quantizeQ81Kernel,
_scaleRowsKernel, _moeWeightedReduceKernel,
_tqRotateQueryKernel, _tqKvAppendKernel, _tqAttentionKernel,
_siluInplaceKernel, _gdnConv1dDecodeKernel, _gdnL2NormPerHeadKernel,
_gdnTileHeadsKernel, _gdnRecurrenceDecodeKernel,
Expand Down Expand Up @@ -3531,6 +3538,8 @@ private void LoadKernelFunctions()
_clearF32Kernel = GetKernelFunc("llm_clear_f32");
_quantizeQ81Kernel = GetKernelFunc("llm_quantize_q8_1");
_bwBaselineKernel = GetKernelFunc("llm_bw_baseline");
_scaleRowsKernel = GetKernelFunc("llm_scale_rows_inplace");
_moeWeightedReduceKernel = GetKernelFunc("llm_moe_weighted_reduce");

// TurboQuant kernels (loaded from the same NVRTC module).
_tqRotateQueryKernel = GetKernelFunc("llm_tq_rotate_query");
Expand Down Expand Up @@ -3755,6 +3764,68 @@ public void AddScaledInPlace(Tensor dst, Tensor src, float scale)
Launch1D(_addScaledKernel, n, args);
}

/// <summary>
/// Issue #129: per-row scalar multiply over a [rows × cols] buffer:
/// <c>buf[i*cols + e] *= scales[i]</c>. The device <paramref name="scales"/> buffer
/// holds one scalar per row. Bit-identical to calling
/// <see cref="ScaleInPlace"/>(row_i, scales[i]) once per row — a single float
/// multiply per element, rounded to float. Used to apply the per-token shared-expert
/// sigmoid gate to the batched shared-expert down output in one launch.
/// </summary>
public void ScaleRowsInPlace(Tensor buf, Tensor scales, int rows, int cols)
{
EnsureImageKernels();
if (!_imageKernelsAvailable)
throw new NotSupportedException("NVRTC is not available; cannot run CUDA image kernels.");
if ((uint)rows > 65535)
throw new ArgumentOutOfRangeException(nameof(rows), rows, "ScaleRowsInPlace: rows must fit the CUDA gridDim.y limit (65535).");

// 2D grid: x walks columns (256-wide blocks), y is the row — the kernel
// recovers (i, e) from block/thread indices with no integer divide.
nint p0 = GetDevPtr(buf);
nint p1 = GetDevPtr(scales);
int p2 = rows, p3 = cols;
nint* args = stackalloc nint[4] { (nint)(&p0), (nint)(&p1), (nint)(&p2), (nint)(&p3) };
uint gridX = (uint)((cols + 255) / 256);
int r = NvrtcInterop.LaunchKernel(_scaleRowsKernel, gridX, (uint)rows, 1, 256, 1, 1, 0, _stream, args, null);
if (r != 0) throw new InvalidOperationException($"cuLaunchKernel(scale_rows) failed: {r}");
}
Comment on lines +3775 to +3792

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

medium

To support the 2D grid mapping and avoid expensive integer division on the GPU, we should launch the kernel with a 2D grid. This also eliminates the potential integer overflow risk from checked((int)total).

    public void ScaleRowsInPlace(Tensor buf, Tensor scales, int rows, int cols)
    {
        EnsureImageKernels();
        if (!_imageKernelsAvailable)
            throw new NotSupportedException("NVRTC is not available; cannot run CUDA image kernels.");

        nint p0 = GetDevPtr(buf);
        nint p1 = GetDevPtr(scales);
        int  p2 = rows, p3 = cols;
        nint* args = stackalloc nint[4] { (nint)(&p0), (nint)(&p1), (nint)(&p2), (nint)(&p3) };
        uint gridX = (uint)((cols + 255) / 256);
        uint gridY = (uint)rows;
        int r = NvrtcInterop.LaunchKernel(_scaleRowsKernel, gridX, gridY, 1, 256, 1, 1, 0, _stream, args, null);
        if (r != 0) throw new InvalidOperationException($"cuLaunchKernel(scale_rows) failed: {r}");
    }


/// <summary>
/// Issue #129: batched MoE top-k weighted scatter-reduce + shared-expert add, in one
/// launch over all N tokens. For each (token i, element e):
/// <c>acc = Σ_k downPartial[(i*na+k)*embDim+e] * weights[i*na+k]; acc += shared[i*embDim+e]; shared[i*embDim+e] = acc;</c>
/// The per-k <c>acc += partial*weight</c> contracts to <c>fmaf</c> (NVRTC fmad=true),
/// one rounding per term, exactly matching the sequential <c>AddScaledInPlace</c> loop;
/// the final <c>acc += shared</c> is a plain add (one rounding), matching the sequential
/// <c>AddInPlace</c>. Routed slots are summed in k=0..na-1 order, shared added last —
/// byte-identical to the per-token <c>Clear + AddScaledInPlace×na + AddInPlace</c>.
/// <paramref name="shared"/> is in/out (must already hold the scaled+rounded shared
/// output per token); each thread owns its (i,e) element so the read-modify-write is
/// race-free.
/// </summary>
public void MoeWeightedReduce(Tensor downPartial, Tensor weights, Tensor shared,
int N, int na, int embDim)
{
EnsureImageKernels();
if (!_imageKernelsAvailable)
throw new NotSupportedException("NVRTC is not available; cannot run CUDA image kernels.");
if ((uint)N > 65535)
throw new ArgumentOutOfRangeException(nameof(N), N, "MoeWeightedReduce: N must fit the CUDA gridDim.y limit (65535).");

// 2D grid: x walks embDim (256-wide blocks), y is the token — the kernel
// recovers (i, e) from block/thread indices with no integer divide/modulo.
nint p0 = GetDevPtr(downPartial);
nint p1 = GetDevPtr(weights);
nint p2 = GetDevPtr(shared);
int p3 = N, p4 = na, p5 = embDim;
nint* args = stackalloc nint[6]
{ (nint)(&p0), (nint)(&p1), (nint)(&p2), (nint)(&p3), (nint)(&p4), (nint)(&p5) };
uint gridX = (uint)((embDim + 255) / 256);
int r = NvrtcInterop.LaunchKernel(_moeWeightedReduceKernel, gridX, (uint)N, 1, 256, 1, 1, 0, _stream, args, null);
if (r != 0) throw new InvalidOperationException($"cuLaunchKernel(moe_weighted_reduce) failed: {r}");
}
Comment on lines +3807 to +3827

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

medium

We should launch MoeWeightedReduce with a 2D grid mapping to avoid expensive integer division and modulo operations on the GPU, while also eliminating the overflow risk of checked((int)total).

    public void MoeWeightedReduce(Tensor downPartial, Tensor weights, Tensor shared,
                                  int N, int na, int embDim)
    {
        EnsureImageKernels();
        if (!_imageKernelsAvailable)
            throw new NotSupportedException("NVRTC is not available; cannot run CUDA image kernels.");

        nint p0 = GetDevPtr(downPartial);
        nint p1 = GetDevPtr(weights);
        nint p2 = GetDevPtr(shared);
        int  p3 = N, p4 = na, p5 = embDim;
        nint* args = stackalloc nint[6]
            { (nint)(&p0), (nint)(&p1), (nint)(&p2), (nint)(&p3), (nint)(&p4), (nint)(&p5) };
        uint gridX = (uint)((embDim + 255) / 256);
        uint gridY = (uint)N;
        int r = NvrtcInterop.LaunchKernel(_moeWeightedReduceKernel, gridX, gridY, 1, 256, 1, 1, 0, _stream, args, null);
        if (r != 0) throw new InvalidOperationException($"cuLaunchKernel(moe_weighted_reduce) failed: {r}");
    }


/// <inheritdoc/>
public void ClampInPlace(Tensor x, float min, float max)
{
Expand Down
50 changes: 50 additions & 0 deletions src/SharpInference.Cuda/CudaKernels.cs
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,56 @@ internal static class CudaKernels
dst[idx] += src[idx] * scale;
}

// ── scale_rows_inplace ────────────────────────────────────────────────────
// Per-row scalar multiply: buf[i*cols + e] *= scales[i]. 2D grid: blockIdx.x
// (× blockDim.x) walks the column e, blockIdx.y is the row i — so there's no
// per-thread integer divide/modulo to recover (i, e). The multiply rounds to
// float exactly like a per-row scale_inplace launch, so the result is
// bit-identical to applying ScaleInPlace(row_i, scales[i]) once per row.
extern ""C"" __global__ void llm_scale_rows_inplace(
float* __restrict__ buf, const float* __restrict__ scales, int rows, int cols)
{
int e = (int)(blockIdx.x * blockDim.x + threadIdx.x);
int i = (int)blockIdx.y;
if (e >= cols || i >= rows) return;
buf[(long)i * cols + e] *= scales[i];
}
Comment on lines +112 to +119

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

medium

Integer division on GPUs is extremely slow. We can completely eliminate it by launching the kernel with a 2D grid mapping, where blockIdx.x maps to the column index and blockIdx.y maps to the row index.

extern "C" __global__ void llm_scale_rows_inplace(
    float* __restrict__ buf, const float* __restrict__ scales, int rows, int cols)
{
    int e = blockIdx.x * blockDim.x + threadIdx.x;
    int i = blockIdx.y;
    if (e >= cols || i >= rows) return;
    buf[(long)i * cols + e] *= scales[i];
}


// ── moe_weighted_reduce ───────────────────────────────────────────────────
// Per-token MoE reduce: for each (token i, element e) sum the na unweighted
// down partials in top-k slot order (k = 0..na-1) with their per-(token,slot)
// weights, then add the already-scaled-and-rounded shared-expert value LAST.
//
// acc = 0
// for k in 0..na-1: acc += downPartial[(i*na+k)*embDim + e] * weights[i*na+k]
// acc += shared[i*embDim + e]
// shared[i*embDim + e] = acc
//
// 2D grid: blockIdx.x (× blockDim.x) walks the element e, blockIdx.y is the
// token i — no per-thread integer divide/modulo to recover (i, e). `shared` is
// in/out — the thread that owns element (i,e) is the only reader and writer, so
// the read-modify-write is race-free. `acc` is a single float register: each
// `acc += p*w` contracts to fmaf under NVRTC's default fmad=true (one rounding
// per term, matching add_scaled_inplace), and the shared add is a plain a+b (one
// rounding, matching add_inplace). Order (routed first, shared last) and per-op
// rounding therefore reproduce the sequential Clear + AddScaledInPlace×na +
// AddInPlace accumulation byte-for-byte.
extern ""C"" __global__ void llm_moe_weighted_reduce(
const float* __restrict__ downPartial, const float* __restrict__ weights,
float* __restrict__ shared, int N, int na, int embDim)
{
int e = (int)(blockIdx.x * blockDim.x + threadIdx.x);
int i = (int)blockIdx.y;
if (e >= embDim || i >= N) return;
float acc = 0.0f;
const float* w = weights + (long)i * na;
const float* p = downPartial + ((long)i * na) * embDim + e;
for (int k = 0; k < na; k++)
acc += p[(long)k * embDim] * w[k];
acc += shared[(long)i * embDim + e];
shared[(long)i * embDim + e] = acc;
}
Comment on lines +140 to +154

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

medium

Similarly, we can optimize llm_moe_weighted_reduce by launching it with a 2D grid mapping to avoid expensive integer division and modulo operations on the GPU.

extern "C" __global__ void llm_moe_weighted_reduce(
    const float* __restrict__ downPartial, const float* __restrict__ weights,
    float* __restrict__ shared, int N, int na, int embDim)
{
    int e = blockIdx.x * blockDim.x + threadIdx.x;
    int i = blockIdx.y;
    if (e >= embDim || i >= N) return;
    float acc = 0.0f;
    const float* w = weights + (long)i * na;
    const float* p = downPartial + ((long)i * na) * embDim + e;
    for (int k = 0; k < na; k++)
        acc += p[(long)k * embDim] * w[k];
    acc += shared[(long)i * embDim + e];
    shared[(long)i * embDim + e] = acc;
}


// ── clamp_inplace ─────────────────────────────────────────────────────────
extern ""C"" __global__ void clamp_inplace(float* __restrict__ x, float lo, float hi, int n)
{
Expand Down
Loading