sundial-core 0.1.0

Linear-programming solver that runs on any GPU as WebGPU compute shaders — restarted PDHG (the PDLP family), with every reported status verified on the CPU in f64.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

struct Params { n: u32, stride: u32, tau: f32, sigma: f32, w: f32, p0: f32, p1: f32, p2: f32 }
@group(0) @binding(0) var<uniform> P: Params;
@group(0) @binding(1) var<storage, read> in_a: array<f32>;
@group(0) @binding(2) var<storage, read> in_b: array<f32>;
@group(0) @binding(3) var<storage, read> in_c: array<f32>;
@group(0) @binding(4) var<storage, read> in_d: array<f32>;
@group(0) @binding(5) var<storage, read> in_e: array<f32>;
@group(0) @binding(6) var<storage, read_write> out_a: array<f32>;
@group(0) @binding(7) var<storage, read_write> out_b: array<f32>;
@group(0) @binding(8) var<storage, read> idx_a: array<u32>;
@group(0) @binding(9) var<storage, read> idx_b: array<u32>;

// x⁺ = clamp(x − τ(c + aᵀy)); x̃ = 2x⁺ − x
// in_a=x  in_b=aty  in_c=c  in_d=lv  in_e=uv  out_a=x_new  out_b=x_tilde
@compute @workgroup_size(256)
fn primal_step(@builtin(global_invocation_id) gid: vec3<u32>) {
    var j = gid.x;
    while (j < P.n) {
        let xj = in_a[j];
        let v = xj - P.tau * (in_c[j] + in_b[j]);
        let xn = clamp(v, in_d[j], in_e[j]);
        out_a[j] = xn;
        out_b[j] = 2.0 * xn - xj;
        j = j + P.stride;
    }
}

// v = y + σ(ax̃); y⁺ = v − σ·clamp(v/σ, lc, uc)
// in_a=y  in_b=axt  in_c=lc  in_d=uc  out_a=y_new
@compute @workgroup_size(256)
fn dual_step(@builtin(global_invocation_id) gid: vec3<u32>) {
    var i = gid.x;
    while (i < P.n) {
        let v = in_a[i] + P.sigma * in_b[i];
        out_a[i] = v - P.sigma * clamp(v / P.sigma, in_c[i], in_d[i]);
        i = i + P.stride;
    }
}

// CSR SpMV, one row per thread, grid-stride over rows.
// idx_a=indptr  idx_b=indices  in_a=values  in_b=x  out_a=out   (P.n = rows)
@compute @workgroup_size(256)
fn spmv(@builtin(global_invocation_id) gid: vec3<u32>) {
    var r = gid.x;
    while (r < P.n) {
        var acc = 0.0;
        for (var k = idx_a[r]; k < idx_a[r + 1u]; k = k + 1u) {
            acc = acc + in_a[k] * in_b[idx_b[k]];
        }
        out_a[r] = acc;
        r = r + P.stride;
    }
}

// out += in (running-sum accumulation for iterate averaging)   in_a=src  out_a=sum
@compute @workgroup_size(256)
fn accum(@builtin(global_invocation_id) gid: vec3<u32>) {
    var i = gid.x;
    while (i < P.n) {
        out_a[i] = out_a[i] + in_a[i];
        i = i + P.stride;
    }
}

// out = in · w (divide running sum by count at check time)   in_a=sum  out_a=avg
@compute @workgroup_size(256)
fn ew_scale(@builtin(global_invocation_id) gid: vec3<u32>) {
    var i = gid.x;
    while (i < P.n) {
        out_a[i] = in_a[i] * P.w;
        i = i + P.stride;
    }
}

// out = a ⊙ b (elementwise multiply; used for row/col unscaling)
// in_a=v  in_b=scale  out_a=out
@compute @workgroup_size(256)
fn ew_mul(@builtin(global_invocation_id) gid: vec3<u32>) {
    var i = gid.x;
    while (i < P.n) {
        out_a[i] = in_a[i] * in_b[i];
        i = i + P.stride;
    }
}