Crate sundial_core 

sundial-core

The Sundial solver. Linear programming that runs on any GPU (Apple, AMD, Intel, NVIDIA) as WebGPU compute shaders via wgpu, natively or in a browser tab. No CUDA.

The algorithm is restarted PDHG, the PDLP family, in f32 WGSL. Operators are matrix-free, so a problem’s constraint matrix never has to exist: a 1,048,576-variable optimal-transport instance solves in about 9.4 seconds on an Apple M4 Pro without materializing a single matrix entry.

The GPU never grades its own homework. Optimal is only returned after an independent CPU-f64 KKT recheck at the exact point being returned. Infeasible and Unbounded are only returned after a CPU-f64 Farkas certificate verifies what the GPU detector nominated. A missed detection is an honest IterationLimit. A false claim is a bug.

use sundial_core::problem::{CsrMatrix, LpProblem, SolveOptions, SolveStatus};

// A carpenter sells tables at $5 and chairs at $3. A table takes 2 planks and
// 1 hour, a chair takes 1 plank and 2 hours, and there are 10 planks and 8
// hours. Sundial minimizes, so maximizing 5t + 3c means minimizing -5t - 3c.
let a = CsrMatrix {
    n_rows: 2,
    n_cols: 2,
    indptr: vec![0, 2, 4],          // where each row starts
    indices: vec![0, 1, 0, 1],      // which column each value sits in
    values: vec![2.0, 1.0, 1.0, 2.0],
};
let lp = LpProblem::new(
    "carpenter".into(),
    a,
    vec![-5.0, -3.0],               // objective, always minimized
    0.0,                            // constant offset
    vec![f64::NEG_INFINITY; 2],     // rows have no lower limit...
    vec![10.0, 8.0],                // ...and these upper limits
    vec![0.0, 0.0],                 // you can't make negative furniture
    vec![f64::INFINITY; 2],         // no ceiling on quantity
)
.unwrap();

let sol = sundial_core::reference::solve(&lp, &SolveOptions::default(), &mut |_| {});
assert_eq!(sol.status, SolveStatus::Optimal);   // 4 tables, 2 chairs, $26

That’s the CPU f64 reference path, which needs no GPU and runs anywhere. For the GPU, swap in gpu::engine::solve_gpu, which is async.

The full worked version, with commentary on every argument:

cargo run --example tiny_lp

Honest limits

Iterate arithmetic is f32, so the headline tolerance is the 1e-4 tier, which is PDLP’s “moderate accuracy.” There’s no presolve. A CPU simplex solver still beats this on small LPs, because the GPU only pays off at scale.

Measured Netlib results, with their caveats attached, are in the repository README.

License

MIT OR Apache-2.0, at your option.

Modules

farkas

Farkas / recession-ray certificate verification. Engines only FLAG candidates (iterate-norm growth streaks at restarts); a status of Infeasible/Unbounded is set exclusively after one of these f64 checks passes — the same certificate honesty as Optimal.

gpu

kkt

linop

CPU-side linear-operator abstraction. PDHG and the KKT certificate only ever need A·x and Aᵀ·y, so a problem may supply kernels instead of a stored matrix (spec: “Matrix-free” problem form).

problem

recover

Integral matching recovery for the taxi demo.

reference

CPU f64 reference implementation of restarted PDHG — the executable statement of this project’s math conventions. The GPU engine mirrors this loop exactly, so this file is the arbiter when the two disagree.

scale

Ruiz equilibration + Pock–Chambolle scaling, per the spec’s conventions: Ā = D_r·A·D_c, c̄ = D_c·c, row bounds ×D_r, col bounds ÷D_c. Unscale: x = D_c·x̄, y = D_r·ȳ.

testgen

Constructed-KKT LP generator: builds (problem, x*, y*) satisfying the KKT conditions exactly, giving an exact optimal oracle with no external solver. Construction: sample sparse A and x*; pick per-column bound status (interior / at-lower / at-upper) and per-row activity (inactive / active upper y*>0 / active lower y*<0); set the active bound equal to (Ax*)ᵢ; choose reduced costs g with signs matching column status; set c = g − Aᵀy*.

transport

Discrete optimal transport between two g×g grids as an LP. Constraint matrix is MATRIX-FREE (row/col sums via TransportOp); the cost VECTOR is materialized, which is what caps the grid size — at g=64 (16.7M variables) the cost vector alone is 67 MiB, so going larger needs the cost computed in-shader. Cells: idx = row*g + col, center ((col+0.5)/g, (row+0.5)/g).

weight

PDLP-style primal-weight balancing. Added because the 1M-variable transport instance stalled under τ = σ, with the dual residual at ~1e-9 while the primal sat at ~1e-4 — badly unbalanced steps, not slow ones.