graphic design tool
git clone https://git.lucas.co/cce-designer.git
feat: Relax's Springs mode runs on the GPU — the first Phase 7 step 4 operator
src/gpu.rs keeps one cce_ui::vk::ComputeDevice per evaluation thread,
opened on first use and kept, under CCE_COMPUTE (auto / cpu / gpu; auto
is CPU under cfg(test) so the suite is the same on every machine).
src/springs.rs is the pattern for every operator after it: one algorithm,
one data layout — flat xyz positions, the rest topology as CSR, pins as
u32 — and two backends held to each other by springs_gpu_matches_cpu
(1.8e-7 worst difference over 1.5k points on an Intel Iris Xe).
Springs is a JACOBI solve now, on both backends. It was Gauss–Seidel
over the edge list in sequence, which no per-point kernel can reproduce;
rather than let the two drift apart, both gather each point's incident
corrections from the pass's starting positions and move once. It
converges about half as fast per iteration, which Iterations already
controls, and the pinned-pull test passes unchanged.
The whole solve is one submission through the API's new run_passes: the
first cut submitted a pass at a time and lost to the CPU at every size
measured, 134k points included, because a round trip costs ~0.5 ms
whatever is inside it. Measured in release, sixteen passes: 1.5k points
cpu 0.25 / gpu 1.5 ms; 15k cpu 2.6 / gpu 3.7 ms; 135k cpu 25 / gpu 14 ms.
GPU_MIN_POINTS = 32k is the auto threshold that follows, and the honest
state of the step: a win for large meshes, a loss for the small ones
most projects have. springs_timing (ignored) is the measurement.
Co-Authored-By: Claude Fable 5.1 <[email protected]>
CLAUDE.md | 52 ++++++++++++
shapeshifter.md | 12 ++-
src/geometry.rs | 40 +++-------
src/gpu.rs | 85 ++++++++++++++++++++
src/main.rs | 140 ++++++++++++++++++++++++++++++++
src/springs.rs | 244 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6 files changed, 543 insertions(+), 30 deletions(-)
diff --git a/CLAUDE.md b/CLAUDE.md
index 5923c79..4b3ed06 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -771,6 +771,58 @@ below native Rust, which is fine for tens of thousands of elements per edit
and wrong for a solver at a million per frame. That is Phase 7's step 4
(WGSL compute through the renderer), not a reason to grow this.
+### GPU compute: the springs solve is the first operator (Phase 7 step 4)
+
+`src/gpu.rs` keeps one `cce_ui::vk::ComputeDevice` per thread, opened on
+first use and kept, so the pipeline cache and the buffers survive from one
+edit to the next; a device costs tens of milliseconds to open and a kernel
+a few to compile, and an operator that paid both per evaluation would lose
+to the CPU every time. `CCE_COMPUTE` decides: unset or `auto` takes the GPU
+when there is one and the operator judges the input big enough; `cpu`
+never opens a device; `gpu` insists, and an operator that cannot get one
+says so through the node-error slot rather than silently taking the CPU
+path. **Under `cfg(test)` auto means CPU**, so the suite is the same on
+every machine and the GPU is exercised only by the tests that ask for it
+by name — the cross-checks. The suite never SETS the variable: libtest
+runs tests in parallel and one that did would race every other test
+reading it (`gpu::parse` is the pure function the choice test covers).
+
+`src/springs.rs` is the pattern every later operator follows: one
+algorithm, one data layout, two backends held to each other by a
+cross-check (`springs_gpu_matches_cpu`, agreement to 1e-4 over 1.5k
+points; it skips with a note where there is no Vulkan). The layout is the
+GPU's — positions as a flat `xyz` array because a `vec3<f32>` in a WGSL
+storage array pads to 16 bytes, the rest topology as CSR with the rest
+length on each incident entry, pins as one `u32` per point — and the CPU
+walks the same arrays in the same order. `solve` chooses the backend; a
+GPU failure in auto mode falls back to the CPU with one stderr note.
+
+**Relax's Springs mode is a JACOBI solve now.** Until 2026-09-24 it was
+Gauss–Seidel over the edge list in sequence, every correction visible to
+the next edge, which no per-point kernel can reproduce; rather than let a
+GPU Jacobi and a CPU Gauss–Seidel drift apart, both run Jacobi: each point
+gathers the corrections of its incident edges from the pass's starting
+positions — half the error toward a free neighbour, all of it toward a
+pinned one — averages them, and moves once. It converges roughly half as
+fast per iteration, which Iterations already controls; the pinned-pull
+test that defines the node's behaviour passes unchanged.
+
+**The whole solve is ONE submission** (`run_passes_over` with a ping-pong
+pair): the topology goes up once, the passes are chained by memory
+barriers with the positions alternating between two device buffers, and
+the result comes back once. The first cut submitted a pass at a time and
+LOST to the CPU at every size measured, 134k points included — a
+submission's round trip is about half a millisecond on an integrated GPU
+whatever the dispatch inside it, and sixteen of them buried a solve that
+takes microseconds. `springs_timing` (ignored; run in release with
+`--ignored --nocapture`) is the measurement, on an Intel Iris Xe, sixteen
+passes: 1.5k points cpu 0.25 ms / gpu 1.5 ms; 15k cpu 2.6 ms / gpu 3.7 ms;
+135k cpu 25 ms / gpu 14 ms, plus ~15 ms of pipeline compile on a device's
+first run. `GPU_MIN_POINTS` (32k) is the auto threshold that follows: the
+GPU is a win for large meshes and a loss for the ones most projects have,
+which is the honest state of step 4 and the reason auto does not simply
+mean GPU.
+
### The volume representation
`src/volume.rs` is a dense signed distance field — `Volume { origin, voxel,
diff --git a/shapeshifter.md b/shapeshifter.md
index cfe4e9d..7b577db 100644
--- a/shapeshifter.md
+++ b/shapeshifter.md
@@ -703,8 +703,16 @@ Nothing else in the workspace loads OpenCL, so the ICD bug leaves with it.
> read the read-write ones back; host-visible mapped buffers, pipelines
> cached by source, every failure an `Err` with naga's diagnostic. Four
> tests run on the machine's Vulkan (Intel Iris Xe here) and skip where
-> there is none. Not yet: a first WGSL operator in this crate held to its
-> CPU twin, and a device kept per evaluation thread.
+> there is none. In this crate: `src/gpu.rs` keeps a device per
+> evaluation thread under `CCE_COMPUTE` (auto / cpu / gpu), and
+> `src/springs.rs` is the first operator — Relax's Springs mode as one
+> Jacobi solve on both backends, held together by `springs_gpu_matches_cpu`.
+> The API grew `run_passes` for it — a whole iterative solve in one
+> submission, ping-ponging on the device — because a pass submitted on its
+> own lost to the CPU at every size; with it the GPU wins from the tens of
+> thousands of points up (135k: 14 ms against 25 ms) and the auto
+> threshold sits at 32k. Next: the rest of the per-point set (Repel,
+> Diffuse, the collision response), which share the pattern.
**Step 4 — GPU compute, through the renderer.** Scripts do not run on the
GPU: no embedded language compiles to GPU code, and none should. Parallel work
diff --git a/src/geometry.rs b/src/geometry.rs
index ec8fbdd..fb66a4a 100644
--- a/src/geometry.rs
+++ b/src/geometry.rs
@@ -1642,8 +1642,10 @@ pub fn resolve_collision_geometry_with_errors(
/// sim state). Vertices carrying `group:<Pin Group>` are pinned: they keep
/// their input position and push everyone else instead — chain a
/// displacement over that group first and this node spreads it through the
-/// surface with a stiffness-shaped falloff. Iterations Gauss–Seidel passes
-/// over the unique edges; Stiffness scales each correction. With no Rest,
+/// surface with a stiffness-shaped falloff. Iterations JACOBI passes over
+/// the rest topology (`crate::springs` — one algorithm on the CPU and the
+/// GPU, the first operator of Phase 7 step 4); Stiffness scales each
+/// correction. With no Rest,
/// an unresolvable Rest, or a Rest whose vertex count differs from the
/// input's, the geometry passes through unchanged — there is nothing
/// coherent to restore toward.
@@ -1718,34 +1720,16 @@ pub fn resolve_relax_geometry_with_errors(
.collect()
};
- let edges: Vec<(usize, usize, f32)> = rest
- .edges()
- .iter()
- .map(|e| {
- let (a, b) = (e[0] as usize, e[1] as usize);
- (a, b, (rest.pos(b) - rest.pos(a)).length())
- })
- .collect();
-
- for _ in 0..iterations {
- for &(a, b, rest_len) in &edges {
- let d = pos[b] - pos[a];
- let len = d.length();
- if len < 1e-6 {
- continue;
- }
- let corr = d * ((len - rest_len) / len * 0.5 * stiffness);
- match (pinned[a], pinned[b]) {
- (false, false) => {
- pos[a] += corr;
- pos[b] -= corr;
- }
- (true, false) => pos[b] -= corr * 2.0,
- (false, true) => pos[a] += corr * 2.0,
- (true, true) => {}
- }
+ let sys = crate::springs::SpringSystem::build(&rest, &pinned);
+ let mut flat: Vec<f32> = pos.iter().flat_map(|v| [v.x, v.y, v.z]).collect();
+ if let Err(e) = crate::springs::solve(&sys, stiffness, iterations, &mut flat) {
+ if ocl_error.is_none() {
+ *ocl_error = Some(format!("{}: {e}", target.name));
}
}
+ for (p, v) in pos.iter_mut().enumerate() {
+ *v = Vec3::new(flat[p * 3], flat[p * 3 + 1], flat[p * 3 + 2]);
+ }
for (p, v) in pos.iter().enumerate() {
geom.set_pos(p, *v);
diff --git a/src/gpu.rs b/src/gpu.rs
new file mode 100644
index 0000000..2115620
--- /dev/null
+++ b/src/gpu.rs
@@ -0,0 +1,85 @@
+//! The evaluation thread's compute device, and the switch that says whether
+//! an operator may use it. Phase 7 step 4 of `shapeshifter.md`.
+//!
+//! One [`ComputeDevice`] per thread, opened on first use and kept, so the
+//! pipeline cache and the buffers survive from one edit to the next: a
+//! device costs tens of milliseconds to open and a kernel a few to compile,
+//! and an operator that paid both per evaluation would lose to the CPU
+//! every time. Thread-local because the device is not `Send`, and because
+//! the app evaluates on its main thread while the thumbnail and export
+//! CLIs evaluate on theirs.
+//!
+//! `CCE_COMPUTE` decides: unset (or `auto`) uses the GPU when one is there
+//! and the operator judges the input large enough to be worth the round
+//! trip; `cpu` never opens a device; `gpu` insists, and an operator that
+//! cannot get one reports it through the node-error slot rather than
+//! silently taking the CPU path. Under `cfg(test)` auto means CPU, so the
+//! suite is deterministic on every machine and the GPU is exercised only by
+//! the tests that ask for it by name — the cross-checks that hold each GPU
+//! operator to its CPU twin.
+
+use cce_ui::vk::ComputeDevice;
+use std::cell::RefCell;
+
+/// What `CCE_COMPUTE` asked for.
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+pub enum Choice {
+ Auto,
+ Cpu,
+ Gpu,
+}
+
+pub fn choice() -> Choice {
+ parse(std::env::var("CCE_COMPUTE").ok().as_deref())
+}
+
+/// `CCE_COMPUTE`'s value to a choice — a pure function of the text, so the
+/// tests can cover it without touching the process environment, which
+/// libtest's parallel tests would race on.
+pub fn parse(value: Option<&str>) -> Choice {
+ let auto = if cfg!(test) { Choice::Cpu } else { Choice::Auto };
+ match value.map(str::trim).map(str::to_ascii_lowercase).as_deref() {
+ Some("cpu") => Choice::Cpu,
+ Some("gpu") => Choice::Gpu,
+ Some("auto") | Some("") | None => auto,
+ Some(other) => {
+ eprintln!("cce-designer: CCE_COMPUTE={other:?} is not cpu, gpu or auto; using auto");
+ auto
+ }
+ }
+}
+
+thread_local! {
+ static DEVICE: RefCell<Option<Result<ComputeDevice, String>>> = const { RefCell::new(None) };
+}
+
+/// Run `f` on this thread's device, opening it on first call. `Err` when
+/// the machine has no usable Vulkan — the same answer every call, since a
+/// failed open is remembered rather than retried per evaluation.
+pub fn with_any_device<R>(f: impl FnOnce(&mut ComputeDevice) -> R) -> Result<R, String> {
+ DEVICE.with(|cell| {
+ let mut slot = cell.borrow_mut();
+ if slot.is_none() {
+ let opened = ComputeDevice::new();
+ match &opened {
+ Ok(d) => eprintln!("cce-designer: compute on {}", d.device_name()),
+ Err(e) => eprintln!("cce-designer: no compute device ({e}); operators run on the CPU"),
+ }
+ *slot = Some(opened);
+ }
+ match slot.as_mut().unwrap() {
+ Ok(d) => Ok(f(d)),
+ Err(e) => Err(e.clone()),
+ }
+ })
+}
+
+/// Whether an operator with `points` elements should take the GPU: the
+/// choice, and in auto mode the operator's own threshold.
+pub fn use_gpu(points: usize, auto_threshold: usize) -> bool {
+ match choice() {
+ Choice::Cpu => false,
+ Choice::Gpu => true,
+ Choice::Auto => points >= auto_threshold,
+ }
+}
diff --git a/src/main.rs b/src/main.rs
index 4035cfb..ffc5786 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -12,6 +12,8 @@ pub mod spatial;
pub mod volume;
pub mod wrangle;
pub mod shapes;
+pub mod gpu;
+pub mod springs;
// Root-level aliases some modules import via `crate::` paths.
#[allow(unused_imports)]
@@ -10599,4 +10601,142 @@ mod tests {
let templates = crate::app::load_fs_tree();
assert!(!templates.children.iter().any(|t| t.node_type == "opencl"), "nodes/opencl.json is gone");
}
+
+ // ---- the springs solve (src/springs.rs): CPU and GPU, one algorithm ----
+
+ /// A rest sphere, its points stretched and stirred, a few pinned: the
+ /// fixture both backends solve.
+ fn springs_fixture(rows: usize, cols: usize) -> (crate::springs::SpringSystem, Vec<f32>, Vec<bool>) {
+ let rest = crate::geometry::sphere_detail(Vec3::ZERO, 0.5, rows, cols);
+ let n = rest.num_points();
+ let pinned: Vec<bool> = (0..n).map(|p| p % 97 == 0).collect();
+ let sys = crate::springs::SpringSystem::build(&rest, &pinned);
+ let pos: Vec<f32> = (0..n)
+ .flat_map(|p| {
+ let x = rest.pos(p);
+ let stretched = x * 1.4 + Vec3::new((p as f32 * 0.37).sin() * 0.05, 0.0, (p as f32 * 0.11).cos() * 0.05);
+ [stretched.x, stretched.y, stretched.z]
+ })
+ .collect();
+ (sys, pos, pinned)
+ }
+
+ /// The CSR is the rest topology twice over — every edge once from each
+ /// end — with the rest length on both entries.
+ #[test]
+ fn springs_system_is_the_rest_topology_in_csr_form() {
+ let rest = crate::geometry::sphere_detail(Vec3::ZERO, 0.5, 6, 8);
+ let sys = crate::springs::SpringSystem::build(&rest, &[]);
+ assert_eq!(sys.n, rest.num_points());
+ assert_eq!(sys.num_edges(), rest.edges().len());
+ assert_eq!(sys.neighbour.len(), 2 * rest.edges().len());
+ for p in 0..sys.n {
+ let (s, e) = (sys.offsets[p] as usize, sys.offsets[p + 1] as usize);
+ assert_eq!(e - s, rest.point_neighbours(p).len(), "point {p}'s incident count is its valence");
+ for i in s..e {
+ let q = sys.neighbour[i] as usize;
+ assert!((sys.rest[i] - (rest.pos(q) - rest.pos(p)).length()).abs() < 1e-6);
+ }
+ }
+ assert!(sys.pinned.iter().all(|&v| v == 0), "no pins asked for, none set");
+ }
+
+ /// Jacobi restores the rest lengths and never moves a pin.
+ #[test]
+ fn springs_cpu_restores_rest_lengths_and_holds_pins() {
+ let (sys, mut pos, pinned) = springs_fixture(12, 16);
+ let before = pos.clone();
+ let error = |pos: &[f32]| -> f32 {
+ let mut worst: f32 = 0.0;
+ for p in 0..sys.n {
+ for e in sys.offsets[p] as usize..sys.offsets[p + 1] as usize {
+ let q = sys.neighbour[e] as usize;
+ let d = Vec3::new(pos[q * 3] - pos[p * 3], pos[q * 3 + 1] - pos[p * 3 + 1], pos[q * 3 + 2] - pos[p * 3 + 2]);
+ worst = worst.max((d.length() - sys.rest[e]).abs() / sys.rest[e]);
+ }
+ }
+ worst
+ };
+ let start = error(&pos);
+ assert!(start > 0.3, "the fixture is stretched: {start}");
+ crate::springs::solve_cpu(&sys, 1.0, 60, &mut pos);
+ let after = error(&pos);
+ assert!(after < start * 0.25, "sixty passes at full stiffness pull the edges toward rest: {start} -> {after}");
+ for p in 0..sys.n {
+ if pinned[p] {
+ assert_eq!(&pos[p * 3..p * 3 + 3], &before[p * 3..p * 3 + 3], "pin {p} moved");
+ }
+ }
+ // Zero stiffness or zero iterations: nothing moves.
+ let mut still = before.clone();
+ crate::springs::solve_cpu(&sys, 0.0, 10, &mut still);
+ assert_eq!(still, before);
+ crate::springs::solve_cpu(&sys, 1.0, 0, &mut still);
+ assert_eq!(still, before);
+ }
+
+ /// The cross-check that holds the GPU to the CPU: the same passes over
+ /// the same arrays agree to floating-point noise. Skips, with a note,
+ /// where the machine has no Vulkan.
+ #[test]
+ fn springs_gpu_matches_cpu() {
+ // Well under the auto threshold: the cross-check asks for the GPU by
+ // name, and a small mesh keeps it quick.
+ let (sys, pos0, _) = springs_fixture(32, 48);
+ let mut cpu = pos0.clone();
+ crate::springs::solve_cpu(&sys, 0.7, 12, &mut cpu);
+ let mut gpu = pos0.clone();
+ let ran = crate::gpu::with_any_device(|dev| crate::springs::solve_gpu(dev, &sys, 0.7, 12, &mut gpu));
+ match ran {
+ Err(e) => {
+ println!("skipping springs_gpu_matches_cpu: {e}");
+ return;
+ }
+ Ok(r) => r.expect("the springs kernel runs"),
+ }
+ let mut worst: f32 = 0.0;
+ for (i, (a, b)) in cpu.iter().zip(&gpu).enumerate() {
+ let d = (a - b).abs();
+ assert!(d < 1e-4, "component {i}: cpu {a} gpu {b}");
+ worst = worst.max(d);
+ }
+ assert!(cpu != pos0, "the solve did something");
+ println!("springs gpu vs cpu: worst component difference {worst:e} over {} points", sys.n);
+ }
+
+ /// `CCE_COMPUTE` parses as written, and under test auto means CPU so the
+ /// suite is the same on every machine.
+ #[test]
+ fn compute_choice_parses_the_variable_and_defaults_to_cpu_under_test() {
+ use crate::gpu::{parse, Choice};
+ assert_eq!(parse(None), Choice::Cpu, "auto is CPU under cfg(test)");
+ assert_eq!(parse(Some("auto")), Choice::Cpu);
+ assert_eq!(parse(Some("gpu")), Choice::Gpu);
+ assert_eq!(parse(Some(" CPU ")), Choice::Cpu, "case-insensitive, trimmed");
+ assert_eq!(parse(Some("banana")), Choice::Cpu, "nonsense is auto, with a note");
+ // The suite never sets the variable: a test that did would race every
+ // other test reading it, and the GPU is exercised by name instead.
+ assert!(std::env::var_os("CCE_COMPUTE").is_none());
+ }
+
+ /// Where the auto threshold sits, and a rough measure of why: the CPU
+ /// and GPU solve timed over a large sphere. Ignored — it is a
+ /// measurement, not an assertion — run with
+ /// `cargo test --release -p cce-designer springs_timing -- --ignored --nocapture`.
+ #[test]
+ #[ignore]
+ fn springs_timing() {
+ for (rows, cols) in [(16, 24), (32, 48), (100, 150), (300, 450)] {
+ let (sys, pos0, _) = springs_fixture(rows, cols);
+ let t = std::time::Instant::now();
+ let mut cpu = pos0.clone();
+ crate::springs::solve_cpu(&sys, 0.7, 16, &mut cpu);
+ let cpu_ms = t.elapsed().as_secs_f64() * 1e3;
+ let mut gpu = pos0.clone();
+ let t = std::time::Instant::now();
+ let ran = crate::gpu::with_any_device(|dev| crate::springs::solve_gpu(dev, &sys, 0.7, 16, &mut gpu));
+ let gpu_ms = t.elapsed().as_secs_f64() * 1e3;
+ println!("{:>7} points, 16 passes: cpu {cpu_ms:8.2} ms gpu {gpu_ms:8.2} ms ({:?})", sys.n, ran.map(|r| r.is_ok()));
+ }
+ }
}
diff --git a/src/springs.rs b/src/springs.rs
new file mode 100644
index 0000000..f29733f
--- /dev/null
+++ b/src/springs.rs
@@ -0,0 +1,244 @@
+//! The edge-length spring solve behind Relax's Springs mode, on the CPU and
+//! on the GPU — the first operator of Phase 7 step 4, and the pattern for
+//! the ones that follow: one algorithm, one data layout, two backends held
+//! to each other by a cross-check.
+//!
+//! It is a JACOBI solve. Every point gathers the corrections of all its
+//! incident edges from the positions as they stood at the start of the
+//! pass, averages them, and moves once; the pass is one dispatch, one
+//! invocation per point. Until 2026-09-24 the CPU solver was Gauss–Seidel
+//! over the edge list in sequence, each correction visible to the next
+//! edge, which no per-point kernel can reproduce — so the parallel form is
+//! the one both backends share, and it is what the CPU runs too, rather
+//! than letting the two drift apart. Jacobi converges more slowly per
+//! iteration (roughly half the rate), which Iterations already controls.
+//!
+//! The layout is the GPU's: positions as a flat `xyz` float array (a
+//! `vec3<f32>` in a WGSL storage array is padded to 16 bytes), the rest
+//! topology as CSR — `offsets[p]..offsets[p + 1]` index the incident edges
+//! of point `p`, each with its neighbour and its rest length — and pins as
+//! one `u32` per point. The CPU solver walks exactly the same arrays in
+//! exactly the same order, so the two agree to floating-point noise.
+
+use crate::detail::Detail;
+use cce_ui::vk::{Binding, ComputeDevice, Kernel};
+use glam::Vec3;
+
+/// The rest topology in CSR form, plus the pins, ready for either backend.
+pub struct SpringSystem {
+ pub n: usize,
+ pub offsets: Vec<u32>,
+ pub neighbour: Vec<u32>,
+ pub rest: Vec<f32>,
+ pub pinned: Vec<u32>,
+}
+
+impl SpringSystem {
+ /// From the REST geometry's unique edges and a per-point pin flag.
+ pub fn build(rest: &Detail, pinned: &[bool]) -> Self {
+ let n = rest.num_points();
+ let mut counts = vec![0u32; n];
+ for e in rest.edges() {
+ counts[e[0] as usize] += 1;
+ counts[e[1] as usize] += 1;
+ }
+ let mut offsets = vec![0u32; n + 1];
+ for p in 0..n {
+ offsets[p + 1] = offsets[p] + counts[p];
+ }
+ let total = offsets[n] as usize;
+ let mut neighbour = vec![0u32; total];
+ let mut rest_len = vec![0f32; total];
+ let mut fill = offsets[..n].to_vec();
+ for e in rest.edges() {
+ let (a, b) = (e[0] as usize, e[1] as usize);
+ let len = (rest.pos(b) - rest.pos(a)).length();
+ let ia = fill[a] as usize;
+ neighbour[ia] = b as u32;
+ rest_len[ia] = len;
+ fill[a] += 1;
+ let ib = fill[b] as usize;
+ neighbour[ib] = a as u32;
+ rest_len[ib] = len;
+ fill[b] += 1;
+ }
+ SpringSystem {
+ n,
+ offsets,
+ neighbour,
+ rest: rest_len,
+ pinned: (0..n).map(|p| u32::from(pinned.get(p).copied().unwrap_or(false))).collect(),
+ }
+ }
+
+ pub fn num_edges(&self) -> usize {
+ self.neighbour.len() / 2
+ }
+}
+
+/// One Jacobi pass on the CPU: `out` from `pos`.
+fn pass_cpu(sys: &SpringSystem, stiffness: f32, pos: &[f32], out: &mut [f32]) {
+ for p in 0..sys.n {
+ let base = p * 3;
+ let x = Vec3::new(pos[base], pos[base + 1], pos[base + 2]);
+ let mut x_new = x;
+ if sys.pinned[p] == 0 {
+ let (start, end) = (sys.offsets[p] as usize, sys.offsets[p + 1] as usize);
+ let mut sum = Vec3::ZERO;
+ for e in start..end {
+ let q = sys.neighbour[e] as usize;
+ let y = Vec3::new(pos[q * 3], pos[q * 3 + 1], pos[q * 3 + 2]);
+ let d = y - x;
+ let len = d.length();
+ if len < 1e-6 {
+ continue;
+ }
+ // Half the error toward a free neighbour (it moves the other
+ // half); all of it toward a pinned one, which does not move.
+ let w = if sys.pinned[q] != 0 { 1.0 } else { 0.5 };
+ sum += d * ((len - sys.rest[e]) / len * w * stiffness);
+ }
+ let count = end - start;
+ if count > 0 {
+ x_new = x + sum / count as f32;
+ }
+ }
+ out[base] = x_new.x;
+ out[base + 1] = x_new.y;
+ out[base + 2] = x_new.z;
+ }
+}
+
+/// `iterations` Jacobi passes on the CPU, in place.
+pub fn solve_cpu(sys: &SpringSystem, stiffness: f32, iterations: usize, pos: &mut Vec<f32>) {
+ let mut out = vec![0f32; pos.len()];
+ for _ in 0..iterations {
+ pass_cpu(sys, stiffness, pos, &mut out);
+ std::mem::swap(pos, &mut out);
+ }
+}
+
+/// The same pass as a WGSL kernel: one invocation per point.
+pub const SPRINGS_WGSL: &str = r#"
+struct Params { stiffness: f32, n: u32, pad0: u32, pad1: u32 }
+@group(0) @binding(0) var<storage, read> pos_in: array<f32>;
+@group(0) @binding(1) var<storage, read_write> pos_out: array<f32>;
+@group(0) @binding(2) var<storage, read> offsets: array<u32>;
+@group(0) @binding(3) var<storage, read> neighbour: array<u32>;
+@group(0) @binding(4) var<storage, read> rest: array<f32>;
+@group(0) @binding(5) var<storage, read> pinned: array<u32>;
+@group(0) @binding(6) var<uniform> params: Params;
+
+@compute @workgroup_size(64)
+fn springs(@builtin(global_invocation_id) id: vec3<u32>) {
+ let p = id.x;
+ if (p >= params.n) { return; }
+ let base = p * 3u;
+ let x = vec3<f32>(pos_in[base], pos_in[base + 1u], pos_in[base + 2u]);
+ var x_new = x;
+ if (pinned[p] == 0u) {
+ let start = offsets[p];
+ let end = offsets[p + 1u];
+ var sum = vec3<f32>(0.0, 0.0, 0.0);
+ for (var e = start; e < end; e = e + 1u) {
+ let q = neighbour[e];
+ let y = vec3<f32>(pos_in[q * 3u], pos_in[q * 3u + 1u], pos_in[q * 3u + 2u]);
+ let d = y - x;
+ let len = length(d);
+ if (len < 1e-6) { continue; }
+ var w = 0.5;
+ if (pinned[q] != 0u) { w = 1.0; }
+ sum = sum + d * ((len - rest[e]) / len * w * params.stiffness);
+ }
+ let count = end - start;
+ if (count > 0u) { x_new = x + sum / f32(count); }
+ }
+ pos_out[base] = x_new.x;
+ pos_out[base + 1u] = x_new.y;
+ pos_out[base + 2u] = x_new.z;
+}"#;
+
+#[repr(C)]
+#[derive(Clone, Copy, bytemuck::Pod, bytemuck::Zeroable)]
+struct Params {
+ stiffness: f32,
+ n: u32,
+ pad0: u32,
+ pad1: u32,
+}
+
+/// `iterations` Jacobi passes on the GPU, in place: ONE submission, the
+/// passes chained by memory barriers and the positions ping-ponging between
+/// two device buffers, so the topology goes up once and the result comes
+/// back once. The first cut submitted a pass at a time — upload, dispatch,
+/// wait, read back, sixteen times — and lost to the CPU at every mesh size
+/// measured, 134k points included: a submission's round trip is about half
+/// a millisecond on an integrated GPU whatever the dispatch inside it, and
+/// the solve itself is far cheaper than that. `springs_timing` has the
+/// numbers for both shapes.
+pub fn solve_gpu(
+ dev: &mut ComputeDevice,
+ sys: &SpringSystem,
+ stiffness: f32,
+ iterations: usize,
+ pos: &mut Vec<f32>,
+) -> Result<(), String> {
+ if sys.n == 0 || iterations == 0 {
+ return Ok(());
+ }
+ let kernel = Kernel::new(SPRINGS_WGSL, "springs");
+ let params = Params { stiffness, n: sys.n as u32, pad0: 0, pad1: 0 };
+ let mut out = vec![0f32; pos.len()];
+ dev.run_passes_over(
+ &kernel,
+ &mut [
+ Binding::input(pos.as_slice()),
+ Binding::rw(out.as_mut_slice()),
+ Binding::input(&sys.offsets),
+ Binding::input(&sys.neighbour),
+ Binding::input(&sys.rest),
+ Binding::input(&sys.pinned),
+ Binding::uniform(¶ms),
+ ],
+ sys.n as u32,
+ iterations as u32,
+ Some((0, 1)),
+ )?;
+ *pos = out;
+ Ok(())
+}
+
+/// Below this many points the CPU is faster; in auto mode it takes them.
+/// From `springs_timing` in release on an Intel Iris Xe, sixteen passes,
+/// one submission: 1.5k points cpu 0.25 ms / gpu 1.5 ms; 15k points cpu
+/// 2.6 ms / gpu 3.7 ms; 135k points cpu 25 ms / gpu 14 ms. Break-even is
+/// in the tens of thousands, and a first run adds ~15 ms of pipeline
+/// compile on top; 32k is where the GPU is clearly ahead on every run
+/// after the first.
+pub const GPU_MIN_POINTS: usize = 32_768;
+
+/// The solve as Relax runs it: the backend `CCE_COMPUTE` and the size
+/// choose, a GPU failure in auto mode falling back to the CPU with a note,
+/// and in forced-GPU mode reported back for the node-error slot.
+pub fn solve(sys: &SpringSystem, stiffness: f32, iterations: usize, pos: &mut Vec<f32>) -> Result<(), String> {
+ if crate::gpu::use_gpu(sys.n, GPU_MIN_POINTS) {
+ let attempt = crate::gpu::with_any_device(|dev| solve_gpu(dev, sys, stiffness, iterations, pos));
+ match attempt {
+ Ok(Ok(())) => return Ok(()),
+ Ok(Err(e)) | Err(e) => {
+ if crate::gpu::choice() == crate::gpu::Choice::Gpu {
+ solve_cpu(sys, stiffness, iterations, pos);
+ return Err(format!("CCE_COMPUTE=gpu but the springs solve could not run there ({e}); solved on the CPU"));
+ }
+ note_fallback_once(&e);
+ }
+ }
+ }
+ solve_cpu(sys, stiffness, iterations, pos);
+ Ok(())
+}
+
+fn note_fallback_once(e: &str) {
+ static ONCE: std::sync::Once = std::sync::Once::new();
+ ONCE.call_once(|| eprintln!("cce-designer: springs solve fell back to the CPU: {e}"));
+}