Deterministic CPU multithreading: supernodal Cholesky, parallel elements, LBM blocks (#50)

* parallel: deterministic thread pool helpers for the solvers (fixed partitions, ordered results)

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* solver: supernodal left-looking Cholesky, multithreaded with a fixed operation order

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* solver: parallel element matrices, stresses, Lanczos block solves and products (sums in the sequential order)

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* explicit: element updates in parallel, nodal forces added in element order (bit-identical to the sequential loop)

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* lbm: collision and streaming by fixed cell blocks in parallel, wall forces in link order (bit-identical)

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* engine: threads setting for the solve commands (session default from MODELCRAFT_THREADS or the cores); bit-identity test on 1, 2 and 7 threads

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* docs: rayon threads in ADR 0001, the parallel crate, supernodal Cholesky

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

---------

Co-authored-by: jyviko <[email protected]>
Co-authored-by: Claude Opus 5.5 <[email protected]>

AGENTS.md

@@ -84,7 +84,8 @@ People trust ModelCraft with their models; a crash loses their work. **This outr
   L7 depends on egui/eframe/winit/wgpu/rfd. **Only
   `crates/geometry` may use `solidcraft-*`.**
 - **Determinism.** The same input gives the same mesh, deck and CPU solver result, bit for bit, on
-  every platform: no hash-order iteration into output, no thread-count-dependent results. GPU
+  every platform: no hash-order iteration into output, no thread-count-dependent results (threads
+  go through `modelcraft-parallel`: fixed blocks, sums in a fixed order, docs/adr/0001). GPU
   solver backends (wgpu, f32) are the one exception: deterministic per device, and checked against
   the CPU reference within a stated tolerance.
 - **Robust geometry.** Meshing decisions (orientation, in-sphere, intersection) use exact or

Cargo.lock

@@ -1901,6 +1901,7 @@ dependencies = [
  "modelcraft-lbm",
  "modelcraft-mesh",
  "modelcraft-model",
+ "modelcraft-parallel",
  "modelcraft-render",
  "modelcraft-solver",
  "serde",
@@ -1915,6 +1916,7 @@ version = "0.1.0"
 dependencies = [
  "modelcraft-fem",
  "modelcraft-model",
+ "modelcraft-parallel",
  "serde",
  "thiserror 2.0.21",
 ]
@@ -1962,6 +1964,7 @@ version = "0.1.0"
 dependencies = [
  "modelcraft-mesh",
  "modelcraft-model",
+ "modelcraft-parallel",
  "thiserror 2.0.21",
 ]
 
@@ -1997,6 +2000,13 @@ dependencies = [
  "thiserror 2.0.21",
 ]
 
+[[package]]
+name = "modelcraft-parallel"
+version = "0.1.0"
+dependencies = [
+ "rayon",
+]
+
 [[package]]
 name = "modelcraft-render"
 version = "0.1.0"
@@ -2013,6 +2023,7 @@ version = "0.1.0"
 dependencies = [
  "modelcraft-fem",
  "modelcraft-model",
+ "modelcraft-parallel",
  "serde",
  "thiserror 2.0.21",
 ]

Cargo.toml

@@ -13,6 +13,7 @@ repository = "https://github.com/jyviko/modelcraft"
 [workspace.dependencies]
 # internal
 modelcraft-model = { path = "crates/model" }
+modelcraft-parallel = { path = "crates/parallel" }
 modelcraft-geometry = { path = "crates/geometry" }
 modelcraft-mesh = { path = "crates/mesh" }
 modelcraft-decks = { path = "crates/decks" }

crates/engine/Cargo.toml

@@ -18,6 +18,7 @@ modelcraft-solver = { workspace = true }
 modelcraft-explicit = { workspace = true }
 modelcraft-fem = { workspace = true }
 modelcraft-lbm = { workspace = true }
+modelcraft-parallel = { workspace = true }
 modelcraft-render = { workspace = true }
 serde = { workspace = true }
 serde_json = { workspace = true }

crates/engine/src/cmd/solve.rs

@@ -5,6 +5,7 @@ use modelcraft_explicit::{ExplicitResult, solve_explicit};
 use modelcraft_lbm::{FluidResult, Obstacle};
 use modelcraft_mesh::surface::{BoundarySurface, SurfaceOptions};
 use modelcraft_model::FluidDomain;
+use modelcraft_parallel as parallel;
 use modelcraft_solver::{
     FrequencyOptions, Limits, Location, MassKind, StaticResult, StepResult as ImplicitResult, solve_frequency, solve_heat, solve_static,
 };
@@ -24,7 +25,7 @@ pub const RUN: CommandSpec = CommandSpec {
 pub const STATIC: CommandSpec = CommandSpec {
     id: "solve.static",
     label: "Solve Static Step",
-    params: r#"{"step": "Load", "max_dofs": 6000000, "max_factor_entries": 400000000} — solves a linear static or coupled temperature–displacement step (the only one when `step` is left out) and keeps the result in the session; the caps may only be lowered"#,
+    params: r#"{"step": "Load", "max_dofs": 6000000, "max_factor_entries": 400000000, "threads": 8} — solves a linear static or coupled temperature–displacement step (the only one when `step` is left out) and keeps the result in the session; the caps may only be lowered; `threads` (1 to 256, by default the session's: `MODELCRAFT_THREADS` or the cores) changes only the time, never the result"#,
     kind: CommandKind::Model,
     run: run_static,
 };
@@ -32,7 +33,7 @@ pub const STATIC: CommandSpec = CommandSpec {
 pub const HEAT: CommandSpec = CommandSpec {
     id: "solve.heat",
     label: "Solve Heat Transfer Step",
-    params: r#"{"step": "Heat", "max_dofs": 6000000, "max_factor_entries": 400000000} — steady or transient heat conduction of a heat transfer step (the only one when `step` is left out): temperatures, heat flux and the heat flow at prescribed temperatures, kept in the session; the caps may only be lowered"#,
+    params: r#"{"step": "Heat", "max_dofs": 6000000, "max_factor_entries": 400000000, "threads": 8} — steady or transient heat conduction of a heat transfer step (the only one when `step` is left out): temperatures, heat flux and the heat flow at prescribed temperatures, kept in the session; the caps may only be lowered; `threads` as for solve.static"#,
     kind: CommandKind::Model,
     run: run_heat,
 };
@@ -40,7 +41,7 @@ pub const HEAT: CommandSpec = CommandSpec {
 pub const FREQUENCY: CommandSpec = CommandSpec {
     id: "solve.frequency",
     label: "Solve Frequency Step",
-    params: r#"{"step": "Modes", "mass": "consistent|lumped", "shift": 0.0, "max_dofs": 6000000, "max_factor_entries": 400000000, "max_iterations": 3000} — natural frequencies and mode shapes of a frequency step (its modes and range); `max_iterations` caps the Lanczos vectors; the result stays in the session"#,
+    params: r#"{"step": "Modes", "mass": "consistent|lumped", "shift": 0.0, "max_dofs": 6000000, "max_factor_entries": 400000000, "max_iterations": 3000, "threads": 8} — natural frequencies and mode shapes of a frequency step (its modes and range); `max_iterations` caps the Lanczos vectors; the result stays in the session; `threads` as for solve.static"#,
     kind: CommandKind::Model,
     run: run_frequency,
 };
@@ -48,7 +49,7 @@ pub const FREQUENCY: CommandSpec = CommandSpec {
 pub const EXPLICIT: CommandSpec = CommandSpec {
     id: "solve.explicit",
     label: "Solve Explicit Dynamic Step",
-    params: r#"{"step": "Impact", "max_dofs": 6000000, "max_increments": 1000000000, "max_output_values": 50000000} — transient response of an explicit dynamic step (the only one when `step` is left out) by central differences; keeps its histories and field frames in the session; the caps may only be lowered (the step's own increment cap applies too)"#,
+    params: r#"{"step": "Impact", "max_dofs": 6000000, "max_increments": 1000000000, "max_output_values": 50000000, "threads": 8} — transient response of an explicit dynamic step (the only one when `step` is left out) by central differences; keeps its histories and field frames in the session; the caps may only be lowered (the step's own increment cap applies too); `threads` as for solve.static"#,
     kind: CommandKind::Model,
     run: run_explicit,
 };
@@ -64,7 +65,7 @@ pub const HISTORY: CommandSpec = CommandSpec {
 pub const FLUID: CommandSpec = CommandSpec {
     id: "solve.fluid",
     label: "Solve Fluid Step",
-    params: r#"{"step": "Flow", "max_cells": 16000000, "max_steps": 10000000, "max_cell_updates": 1e12, "max_frame_bytes": 2e9} — runs a fluid step with the lattice Boltzmann solver (the only fluid step when `step` is left out): obstacles are triangulated from the open design or taken from imported surfaces; returns the time step, relaxation time, Mach number, mass and the final force on every wall; the caps may only be lowered"#,
+    params: r#"{"step": "Flow", "max_cells": 16000000, "max_steps": 10000000, "max_cell_updates": 1e12, "max_frame_bytes": 2e9, "threads": 8} — runs a fluid step with the lattice Boltzmann solver (the only fluid step when `step` is left out): obstacles are triangulated from the open design or taken from imported surfaces; returns the time step, relaxation time, Mach number, mass and the final force on every wall; the caps may only be lowered; `threads` as for solve.static"#,
     kind: CommandKind::Model,
     run: run_fluid,
 };
@@ -99,6 +100,18 @@ fn cap<T: Copy + TryInto<u64> + TryFrom<u64> + std::fmt::Display>(p: &Value, key
     }
 }
 
+/// The `threads` parameter: 1 to [`parallel::MAX_THREADS`], by default the session's.
+fn threads(s: &Session, p: &Value) -> Result<usize> {
+    match p.get("threads") {
+        None => Ok(s.limits.threads.clamp(1, parallel::MAX_THREADS)),
+        Some(v) => v
+            .as_u64()
+            .and_then(|n| usize::try_from(n).ok())
+            .filter(|n| (1..=parallel::MAX_THREADS).contains(n))
+            .ok_or_else(|| EngineError::Params(format!("`threads` must be a whole number from 1 to {}", parallel::MAX_THREADS))),
+    }
+}
+
 fn limits(p: &Value) -> Result<Limits> {
     Ok(Limits {
         max_dofs: cap(p, "max_dofs", Limits::MAX_DOFS)?,
@@ -155,7 +168,8 @@ fn smallest(values: impl Iterator<Item = (u32, f64)>) -> Value {
 fn run_static(s: &mut Session, p: &Value) -> Result<Value> {
     let limits = limits(p)?;
     let step = opt::<String>(p, "step")?;
-    let r = solve_static(&s.model, step.as_deref(), &limits)?;
+    let model = &s.model;
+    let r = parallel::with_threads(threads(s, p)?, || solve_static(model, step.as_deref(), &limits))?;
 
     let summary = json!({
         "step": r.step,
@@ -182,7 +196,8 @@ fn run_static(s: &mut Session, p: &Value) -> Result<Value> {
 fn run_heat(s: &mut Session, p: &Value) -> Result<Value> {
     let limits = limits(p)?;
     let step = opt::<String>(p, "step")?;
-    let r = solve_heat(&s.model, step.as_deref(), &limits)?;
+    let model = &s.model;
+    let r = parallel::with_threads(threads(s, p)?, || solve_heat(model, step.as_deref(), &limits))?;
     let t = r.temperature();
     let summary = json!({
         "step": r.step,
@@ -206,7 +221,8 @@ fn run_frequency(s: &mut Session, p: &Value) -> Result<Value> {
     let step = opt::<String>(p, "step")?;
     let mass = opt::<MassKind>(p, "mass")?.unwrap_or_default();
     let shift = opt::<f64>(p, "shift")?;
-    let r = solve_frequency(&s.model, step.as_deref(), &FrequencyOptions { mass, shift }, &limits)?;
+    let model = &s.model;
+    let r = parallel::with_threads(threads(s, p)?, || solve_frequency(model, step.as_deref(), &FrequencyOptions { mass, shift }, &limits))?;
     let summary = json!({
         "step": r.step,
         "dofs": r.stats.dofs,
@@ -233,7 +249,8 @@ fn run_explicit(s: &mut Session, p: &Value) -> Result<Value> {
         max_output_values: cap(p, "max_output_values", ExplicitLimits::MAX_OUTPUT_VALUES)?,
     };
     let step = opt::<String>(p, "step")?;
-    let r = solve_explicit(&s.model, step.as_deref(), &limits)?;
+    let model = &s.model;
+    let r = parallel::with_threads(threads(s, p)?, || solve_explicit(model, step.as_deref(), &limits))?;
     let summary = explicit_summary(&r);
     s.store_result(StepResult::Explicit(r));
     Ok(summary)
@@ -333,7 +350,8 @@ fn run_fluid(s: &mut Session, p: &Value) -> Result<Value> {
     let step = opt::<String>(p, "step")?;
     let (_, domain) = modelcraft_lbm::fluid_step(&s.model, step.as_deref())?;
     let obstacles = obstacles(s, domain)?;
-    let r = modelcraft_lbm::solve_fluid(&s.model, step.as_deref(), &obstacles, &limits)?;
+    let model = &s.model;
+    let r = parallel::with_threads(threads(s, p)?, || modelcraft_lbm::solve_fluid(model, step.as_deref(), &obstacles, &limits))?;
     let summary = fluid_summary(&r);
     s.store_result(StepResult::Fluid(r));
     Ok(summary)
@@ -756,6 +774,13 @@ pub(crate) mod tests {
             json!({"step": "Modes", "nodes": [1], "mode": 1e300}),
             json!({"step": "Bend", "nodes": [0, 99999999, 4294967296u64]}),
             json!({"step": "Bend", "nodes": vec![1; 20_000]}),
+            json!({"step": "Bend", "threads": 0}),
+            json!({"step": "Bend", "threads": -3}),
+            json!({"step": "Bend", "threads": 257}),
+            json!({"step": "Bend", "threads": 1e300}),
+            json!({"step": "Bend", "threads": "all"}),
+            json!({"step": "Modes", "threads": 18446744073709551615u64}),
+            json!({"step": "Modes", "threads": 2.5}),
         ];
         for id in ["solve.static", "solve.frequency", "solve.probe"] {
             for p in &hostile {
@@ -896,6 +921,9 @@ pub(crate) mod tests {
             json!({"step": "Shear", "forces": vec!["y_min"; 1000]}),
             json!({"step": "Shear", "forces": "y_min"}),
             json!({"step": "Shear", "nodes": [1]}),
+            json!({"step": "Shear", "threads": 0}),
+            json!({"step": "Shear", "threads": 100000}),
+            json!({"step": "Shear", "threads": null}),
         ];
         for id in ["solve.fluid", "solve.probe", "solve.static", "solve.frequency"] {
             for p in &hostile {
@@ -972,6 +1000,69 @@ pub(crate) mod tests {
         s.run_script(&script).unwrap();
     }
 
+    /// The stored result of a step, hashed bit for bit (Debug prints every float exactly).
+    fn result_hash(s: &Session, step: &str) -> u64 {
+        use std::hash::{Hash, Hasher};
+        let mut h = std::hash::DefaultHasher::new();
+        format!("{:?}", s.result(step).unwrap()).hash(&mut h);
+        h.finish()
+    }
+
+    /// Implicit (static, frequency), explicit and lattice Boltzmann solves large enough for every
+    /// parallel path (factorization tasks and tiles, element rounds, cell blocks) give the same
+    /// bits on 1, 2 and 7 threads.
+    #[test]
+    fn results_do_not_depend_on_the_thread_count() {
+        let solve_all = |threads: usize| -> Vec<u64> {
+            let mut s = Session::default();
+            explicit_cantilever(&mut s, 10.0, [30, 6]);
+            let script = json!([
+                {"command": "model.load.create", "params": {"name": "push", "type": "force", "target": {"set": "x30"}, "vector": [0, 0.3, -1]}},
+                {"command": "model.step.create", "params": {"name": "Bend", "procedure": {"type": "static"}, "boundary_conditions": ["clamp"], "loads": ["push"]}},
+            ]);
+            s.run_script(&script).unwrap();
+            s.execute("solve.static", &json!({"step": "Bend", "threads": threads})).unwrap();
+            s.execute("solve.frequency", &json!({"step": "Modes", "threads": threads})).unwrap();
+            let mut hashes = vec![result_hash(&s, "Bend"), result_hash(&s, "Modes")];
+            let step = json!({
+                "type": "explicit_dynamic", "end_time": 0.15, "output_interval": 0.05,
+                "initial_velocities": [{"set": "x30", "velocity": [0, 0, -2]}], "hourglass": {"type": "stiffness"}
+            });
+            let script = json!([
+                {"command": "model.section.edit", "params": {"name": "beam", "set": "all", "material": "Mat", "kind": {"type": "solid"}, "formulation": {"integration": "reduced"}}},
+                {"command": "model.step.create", "params": {"name": "Hit", "procedure": step, "boundary_conditions": ["clamp"]}},
+            ]);
+            s.run_script(&script).unwrap();
+            s.execute("solve.explicit", &json!({"step": "Hit", "threads": threads})).unwrap();
+            hashes.push(result_hash(&s, "Hit"));
+            let mut f = Session::default();
+            let periodic = json!({"type": "periodic"});
+            let script = json!([
+                {"command": "model.material.create", "params": {"name": "Oil", "density": 1.0, "fluid": {"type": "newtonian", "viscosity": 0.1}}},
+                {"command": "model.fluid_domain.create", "params": {"name": "Gap", "material": "Oil", "min": [0, 0, 0], "max": [5, 1, 2], "spacing": 0.125,
+                    "faces": {"z_min": periodic, "z_max": periodic}}},
+                {"command": "model.fluid_domain.boundary", "params": {"domain": "Gap", "at": "y_max", "condition": {"type": "moving_wall", "velocity": [0.1, 0, 0.02]}}},
+                {"command": "model.step.create", "params": {"name": "Shear", "procedure": {"type": "fluid", "domain": "Gap", "end_time": 2, "frames": 2}}},
+            ]);
+            f.run_script(&script).unwrap();
+            f.execute("solve.fluid", &json!({"step": "Shear", "threads": threads})).unwrap();
+            hashes.push(result_hash(&f, "Shear"));
+            hashes
+        };
+        let one = solve_all(1);
+        for threads in [2, 7] {
+            assert_eq!(solve_all(threads), one, "{threads} threads");
+        }
+        let mut s = Session::default();
+        cantilever(&mut s);
+        for bad in [json!(0), json!(257), json!(-1), json!(1.5), json!("8")] {
+            let e = s.execute("solve.static", &json!({"threads": bad})).unwrap_err().to_string();
+            assert!(e.contains("`threads` must be a whole number from 1 to 256"), "{e}");
+        }
+        s.limits.threads = 3;
+        s.execute("solve.static", &json!({})).unwrap();
+    }
+
     /// Upward zero crossings of a sampled signal, linearly interpolated.
     fn crossings(time: &[f64], x: &[f64]) -> Vec<f64> {
         let mut out = Vec::new();

crates/engine/src/lib.rs

@@ -135,12 +135,16 @@ pub struct Limits {
     pub max_steiner: usize,
     /// `mesh.tet.generate`'s `max_tets`.
     pub max_tets: usize,
+    /// Threads the solvers use (the `solve.*` commands' `threads`): by default
+    /// `MODELCRAFT_THREADS` or the available cores, 1 to [`modelcraft_parallel::MAX_THREADS`].
+    /// Results do not depend on it.
+    pub threads: usize,
 }
 
 impl Default for Limits {
     fn default() -> Self {
         let cdt = modelcraft_mesh::cdt::CdtOptions::default();
-        Limits { max_steiner: cdt.max_steiner, max_tets: cdt.max_tets }
+        Limits { max_steiner: cdt.max_steiner, max_tets: cdt.max_tets, threads: modelcraft_parallel::default_threads() }
     }
 }
 
@@ -469,7 +473,7 @@ mod tests {
         // Small work caps, so input that would mesh up to the production caps (needle triangles
         // need Steiner points all along them) is refused at once. The caps themselves are what
         // `tet_generate_honours_the_session_caps` checks.
-        let capped = || Session { limits: Limits { max_steiner: 64, max_tets: 100_000 }, ..Session::default() };
+        let capped = || Session { limits: Limits { max_steiner: 64, max_tets: 100_000, ..Limits::default() }, ..Session::default() };
         for spec in commands() {
             for p in &hostile {
                 let mut s = capped();

crates/explicit/Cargo.toml

@@ -11,5 +11,6 @@ workspace = true
 [dependencies]
 modelcraft-model = { workspace = true }
 modelcraft-fem = { workspace = true }
+modelcraft-parallel = { workspace = true }
 thiserror = { workspace = true }
 serde = { workspace = true }

crates/explicit/src/elements.rs

@@ -372,20 +372,111 @@ pub fn hex_scales(
     HexScales { size, coefficient, omega2 }
 }
 
-/// Evaluate the internal forces at the motion `m` into `ws`, updating the element states, and
-/// return the smallest critical increment of the elements still in the analysis.
-pub fn internal_forces(sys: &System, m: &Motion, state: &mut State, ws: &mut Workspace) -> Result<f64> {
-    for f in [&mut ws.stress_force, &mut ws.hourglass_force, &mut ws.viscous_force] {
-        f.fill(0.0);
+/// Elements per parallel task.
+const BLOCK: usize = 128;
+
+/// Tasks whose element forces are computed before they are added (bounds the memory of the
+/// element forces held at once).
+const ROUND: usize = 256;
+
+/// What one element gives an evaluation.
+enum Outcome<const N: usize> {
+    /// Eroded earlier: nothing.
+    Skipped,
+    /// Eroded now.
+    Eroded(Failure),
+    Forces(Forces<N>),
 }
-    state.eroded.clear();
-    let mut critical = f64::INFINITY;
-    let hourglass = modelcraft_model::Hourglass { kind: sys.hourglass, coefficient: sys.hourglass_coefficient };
-    for h in &sys.hexes {
-        let (Some(mat), Some(es)) = (sys.materials.get(h.material), state.elements.get_mut(h.index)) else { continue };
+
+/// An element's nodal forces by origin and its critical increment.
+struct Forces<const N: usize> {
+    hourglass: Option<[[f64; 3]; N]>,
+    stress: [[f64; 3]; N],
+    viscous: Option<[[f64; 3]; N]>,
+    critical: f64,
+}
+
+/// A block of consecutive elements of one kind (positions `items` in their list) with the
+/// slice of element states that holds theirs, from state index `base` on.
+struct Piece<'s, T> {
+    items: std::ops::Range<usize>,
+    base: usize,
+    states: &'s mut [ElementState],
+    out: Vec<Result<T>>,
+}
+
+/// Cut the element states into one slice per block of [`BLOCK`] elements of a kind whose state
+/// indices (`index(i)` for item i) increase, as [`System`] builds them.
+fn pieces<'s, T>(states: &'s mut [ElementState], count: usize, index: &impl Fn(usize) -> usize) -> Result<Vec<Piece<'s, T>>> {
+    let ranges: Vec<std::ops::Range<usize>> = (0..count.div_ceil(BLOCK)).map(|k| k * BLOCK..((k + 1) * BLOCK).min(count)).collect();
+    let mut out = Vec::with_capacity(ranges.len());
+    let mut rest = states;
+    let mut at = 0;
+    for (k, items) in ranges.iter().enumerate() {
+        let base = index(items.start);
+        let end = ranges.get(k + 1).map_or(at + rest.len(), |r| index(r.start));
+        if base < at || end < base || end - at > rest.len() {
+            return Err(Error::Internal("element states out of order".into()));
+        }
+        let (_, tail) = std::mem::take(&mut rest).split_at_mut(base - at);
+        let (mine, tail) = tail.split_at_mut(end - base);
+        out.push(Piece { items: items.clone(), base, states: mine, out: Vec::new() });
+        rest = tail;
+        at = end;
+    }
+    Ok(out)
+}
+
+/// `f(i, state)` for every element i of a kind, in parallel by blocks and a round of blocks at a
+/// time; the outcomes go to `apply` in element order, so whatever it sums is summed as a
+/// sequential loop would. The first error in element order stops the evaluation.
+fn run<T: Send>(
+    states: &mut [ElementState],
+    count: usize,
+    index: impl Fn(usize) -> usize + Sync,
+    f: impl Fn(usize, &mut ElementState) -> Result<T> + Sync + Send,
+    mut apply: impl FnMut(usize, T) -> Result<()>,
+) -> Result<()> {
+    let mut all = pieces::<T>(states, count, &index)?;
+    for round in all.chunks_mut(ROUND) {
+        modelcraft_parallel::for_each_mut(round, |_, piece| {
+            let Piece { items, base, states, out } = piece;
+            *out = items
+                .clone()
+                .map(|i| match states.get_mut(index(i).wrapping_sub(*base)) {
+                    Some(es) => f(i, es),
+                    None => Err(Error::Internal("element state".into())),
+                })
+                .collect();
+        });
+        for piece in round.iter_mut() {
+            for (i, outcome) in piece.items.clone().zip(std::mem::take(&mut piece.out)) {
+                apply(i, outcome?)?;
+            }
+        }
+    }
+    Ok(())
+}
+
+/// Add an element's forces into the workspace and return its critical increment.
+fn add<const N: usize>(ws: &mut Workspace, nodes: &[usize; N], f: &Forces<N>) -> f64 {
+    if let Some(h) = &f.hourglass {
+        scatter(&mut ws.hourglass_force, nodes, h);
+    }
+    scatter(&mut ws.stress_force, nodes, &f.stress);
+    if let Some(q) = &f.viscous {
+        scatter(&mut ws.viscous_force, nodes, q);
+    }
+    f.critical
+}
+
+/// One hexahedron's update and forces.
+fn hexahedron(sys: &System, h: &crate::system::Hex, m: &Motion, es: &mut ElementState) -> Result<Outcome<8>> {
+    let Some(mat) = sys.materials.get(h.material) else { return Ok(Outcome::Skipped) };
     if es.failed.is_some() {
-            continue;
+        return Ok(Outcome::Skipped);
     }
+    let hourglass = modelcraft_model::Hourglass { kind: sys.hourglass, coefficient: sys.hourglass_coefficient };
     let e = Nodal::gather(sys, &h.nodes, m);
     let x_now = e.positions(0.0);
     let (mid, now) = if sys.nonlinear {
@@ -401,8 +492,7 @@ pub fn internal_forces(sys: &System, m: &Motion, state: &mut State, ws: &mut Wor
     let (rho, c) = current(mat, h.kinematics.volume, now.volume);
     update(mat, es, &q, &d, d[0] + d[1] + d[2], m.dt, rho, mid.volume);
     if let Some(failure) = es.failed {
-            state.eroded.push((h.index, failure));
-            continue;
+        return Ok(Outcome::Eroded(failure));
     }
     // Hourglass forces with the current shape vectors and size (fixed without nonlinear
     // geometry).
@@ -433,19 +523,16 @@ pub fn internal_forces(sys: &System, m: &Motion, state: &mut State, ws: &mut Wor
             }
         }
     }
-        scatter(&mut ws.hourglass_force, &h.nodes, &f);
-        let element = Element { now: &now, mid: &mid, nodes: &h.nodes, density: rho, wave_speed: c, length: 2.0 * c / scales.omega2.sqrt() };
-        critical = critical.min(element.finish(ws, &es.point.stress, &e.v, m.bulk));
+    let element = Element { now: &now, mid: &mid, density: rho, wave_speed: c, length: 2.0 * c / scales.omega2.sqrt() };
+    let mut forces = element.finish(&es.point.stress, &e.v, m.bulk);
+    forces.hourglass = Some(f);
+    Ok(Outcome::Forces(forces))
 }
-    if !sys.tets.is_empty() {
-        // Pass 1: each tetrahedron's increment, and each node's mean volumetric increment.
-        ws.theta.fill(0.0);
-        ws.tet_volume.fill(0.0);
-        for (k, t) in sys.tets.iter().enumerate() {
-            let Some(slot) = ws.tet_increment.get_mut(k) else { continue };
-            *slot = None;
+
+/// A tetrahedron's first pass: its increment (none when it has eroded).
+fn tetrahedron_increment(sys: &System, t: &crate::system::Tet, m: &Motion, state: &State) -> Result<Option<TetIncrement>> {
     if state.elements.get(t.index).is_none_or(|e| e.failed.is_some()) {
-                continue;
+        return Ok(None);
     }
     let e = Nodal::gather(sys, &t.nodes, m);
     let (mid, now) = if sys.nonlinear {
@@ -455,27 +542,20 @@ pub fn internal_forces(sys: &System, m: &Motion, state: &mut State, ws: &mut Wor
         (t.kinematics, t.kinematics)
     };
     let g = mid.gradient(&e.du);
-            for node in t.nodes {
-                if let (Some(x), Some(w)) = (ws.theta.get_mut(node), ws.tet_volume.get_mut(node)) {
-                    *x += mid.volume * (g[0][0] + g[1][1] + g[2][2]);
-                    *w += mid.volume;
+    Ok(Some(TetIncrement { g, mid, now, v: e.v }))
 }
-            }
-            *slot = Some(TetIncrement { g, mid, now, v: e.v });
-        }
-        for (x, w) in ws.theta.iter_mut().zip(&ws.tet_volume) {
-            if *w > 0.0 {
-                *x /= w;
-            }
-        }
-        // Pass 2: each tetrahedron's update with the mean of its nodes' volumetric increments.
-        for (k, t) in sys.tets.iter().enumerate() {
-            let (Some(mat), Some(es), Some(Some(inc))) =
-                (sys.materials.get(t.material), state.elements.get_mut(t.index), ws.tet_increment.get(k).copied())
-            else {
-                continue;
-            };
-            let averaged = t.nodes.iter().map(|n| ws.theta.get(*n).copied().unwrap_or(0.0)).sum::<f64>() / 4.0;
+
+/// A tetrahedron's second pass: its update with the mean of its nodes' volumetric increments.
+fn tetrahedron(
+    sys: &System,
+    t: &crate::system::Tet,
+    inc: Option<TetIncrement>,
+    theta: &[f64],
+    m: &Motion,
+    es: &mut ElementState,
+) -> Result<Outcome<4>> {
+    let (Some(mat), Some(inc)) = (sys.materials.get(t.material), inc) else { return Ok(Outcome::Skipped) };
+    let averaged = t.nodes.iter().map(|n| theta.get(*n).copied().unwrap_or(0.0)).sum::<f64>() / 4.0;
     let q = rotation_of(sys, mat, es, &inc.g, t.label, m.time)?;
     let (d, _) = rotation::split(&inc.g);
     let thermal = t.thermal_rate * m.dt;
@@ -483,12 +563,97 @@ pub fn internal_forces(sys: &System, m: &Motion, state: &mut State, ws: &mut Wor
     let (rho, c) = current(mat, t.kinematics.volume, inc.now.volume);
     update(mat, es, &q, &d, averaged - 3.0 * thermal, m.dt, rho, inc.mid.volume);
     if let Some(failure) = es.failed {
-                state.eroded.push((t.index, failure));
-                continue;
+        return Ok(Outcome::Eroded(failure));
     }
     let omega2 = if sys.nonlinear { frequency_bound(&inc.now, inc.now.gradient_norm(), mat.modulus, &t.mass) } else { t.omega2 };
-            let element = Element { now: &inc.now, mid: &inc.mid, nodes: &t.nodes, density: rho, wave_speed: c, length: 2.0 * c / omega2.sqrt() };
-            critical = critical.min(element.finish(ws, &es.point.stress, &inc.v, m.bulk));
+    let element = Element { now: &inc.now, mid: &inc.mid, density: rho, wave_speed: c, length: 2.0 * c / omega2.sqrt() };
+    Ok(Outcome::Forces(element.finish(&es.point.stress, &inc.v, m.bulk)))
+}
+
+/// Evaluate the internal forces at the motion `m` into `ws`, updating the element states, and
+/// return the smallest critical increment of the elements still in the analysis.
+///
+/// The elements are updated in parallel (each owns its state); their nodal forces are added in
+/// element order, hexahedra first, so every sum is the one of a sequential loop and the result
+/// does not depend on the number of threads.
+pub fn internal_forces(sys: &System, m: &Motion, state: &mut State, ws: &mut Workspace) -> Result<f64> {
+    for f in [&mut ws.stress_force, &mut ws.hourglass_force, &mut ws.viscous_force] {
+        f.fill(0.0);
+    }
+    state.eroded.clear();
+    let mut critical = f64::INFINITY;
+    let State { elements, eroded } = state;
+    let hexes = &sys.hexes;
+    run(
+        elements,
+        hexes.len(),
+        |i| hexes.get(i).map_or(usize::MAX, |h| h.index),
+        |i, es| match hexes.get(i) {
+            Some(h) => hexahedron(sys, h, m, es),
+            None => Ok(Outcome::Skipped),
+        },
+        |i, outcome| {
+            let Some(h) = hexes.get(i) else { return Ok(()) };
+            match outcome {
+                Outcome::Skipped => {}
+                Outcome::Eroded(failure) => eroded.push((h.index, failure)),
+                Outcome::Forces(f) => critical = critical.min(add(ws, &h.nodes, &f)),
+            }
+            Ok(())
+        },
+    )?;
+    if !sys.tets.is_empty() {
+        let tets = &sys.tets;
+        // Pass 1: each tetrahedron's increment, and each node's mean volumetric increment.
+        let increments = {
+            let st: &State = state;
+            modelcraft_parallel::map_grained(tets.len(), BLOCK, |k| match tets.get(k) {
+                Some(t) => tetrahedron_increment(sys, t, m, st),
+                None => Ok(None),
+            })
+        };
+        ws.theta.fill(0.0);
+        ws.tet_volume.fill(0.0);
+        for ((t, slot), inc) in tets.iter().zip(ws.tet_increment.iter_mut()).zip(increments) {
+            *slot = inc?;
+            let Some(inc) = slot else { continue };
+            for node in t.nodes {
+                if let (Some(x), Some(w)) = (ws.theta.get_mut(node), ws.tet_volume.get_mut(node)) {
+                    *x += inc.mid.volume * (inc.g[0][0] + inc.g[1][1] + inc.g[2][2]);
+                    *w += inc.mid.volume;
+                }
+            }
+        }
+        for (x, w) in ws.theta.iter_mut().zip(&ws.tet_volume) {
+            if *w > 0.0 {
+                *x /= w;
+            }
+        }
+        // Pass 2: each tetrahedron's update.
+        let State { elements, eroded } = state;
+        let (theta, increments) = (&ws.theta, &ws.tet_increment);
+        let mut forces: Vec<(usize, Forces<4>)> = Vec::new();
+        run(
+            elements,
+            tets.len(),
+            |k| tets.get(k).map_or(usize::MAX, |t| t.index),
+            |k, es| match tets.get(k) {
+                Some(t) => tetrahedron(sys, t, increments.get(k).copied().flatten(), theta, m, es),
+                None => Ok(Outcome::Skipped),
+            },
+            |k, outcome| {
+                match (tets.get(k), outcome) {
+                    (Some(t), Outcome::Eroded(failure)) => eroded.push((t.index, failure)),
+                    (Some(_), Outcome::Forces(f)) => forces.push((k, f)),
+                    _ => {}
+                }
+                Ok(())
+            },
+        )?;
+        for (k, f) in &forces {
+            if let Some(t) = tets.get(*k) {
+                critical = critical.min(add(ws, &t.nodes, f));
+            }
         }
     }
     for r in &sys.rods {
@@ -524,21 +689,21 @@ struct TetIncrement {
 struct Element<'a, const N: usize> {
     now: &'a UniformStrain<N>,
     mid: &'a UniformStrain<N>,
-    nodes: &'a [usize; N],
     density: f64,
     wave_speed: f64,
     length: f64,
 }
 
 impl<const N: usize> Element<'_, N> {
-    /// Add the stress and bulk viscosity forces and return the critical increment.
-    fn finish(&self, ws: &mut Workspace, sigma: &[f64; 6], v: &[[f64; 3]; N], bulk: &BulkViscosity) -> f64 {
+    /// The stress and bulk viscosity forces and the critical increment.
+    fn finish(&self, sigma: &[f64; 6], v: &[[f64; 3]; N], bulk: &BulkViscosity) -> Forces<N> {
         let rate = self.mid.divergence(v);
         let q = bulk_pressure(rate, self.density, self.length, self.wave_speed, bulk);
-        scatter(&mut ws.stress_force, self.nodes, &self.now.forces(sigma));
-        if q != 0.0 {
-            scatter(&mut ws.viscous_force, self.nodes, &self.now.pressure_forces(q));
+        Forces {
+            hourglass: None,
+            stress: self.now.forces(sigma),
+            viscous: (q != 0.0).then(|| self.now.pressure_forces(q)),
+            critical: critical_increment(rate, self.length, self.wave_speed, bulk),
         }
-        critical_increment(rate, self.length, self.wave_speed, bulk)
     }
 }

crates/explicit/src/lib.rs

@@ -24,9 +24,10 @@
 //! docs/solvers/literature.md, Parts F and J). Without it the gradients stay those of the
 //! undeformed mesh (small displacements).
 //!
-//! Everything is sequential, with elements in label order and no parallel reductions, so the
-//! same model gives the same numbers bit for bit; sizes, increments and output are capped by
-//! [`Limits`].
+//! Elements are updated on threads ([`modelcraft_parallel`]), each owning its state, and their
+//! nodal forces are added in label order as a sequential loop adds them, with no parallel
+//! reductions, so the same model gives the same numbers bit for bit on any number of threads;
+//! sizes, increments and output are capped by [`Limits`].
 #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)]
 #![forbid(unsafe_code)]
 // Small fixed-size vector kernels read as the textbook formulas with index loops.

crates/lbm/Cargo.toml

@@ -11,4 +11,5 @@ workspace = true
 [dependencies]
 modelcraft-model = { workspace = true }
 modelcraft-mesh = { workspace = true }
+modelcraft-parallel = { workspace = true }
 thiserror = { workspace = true }

crates/lbm/src/flow.rs

@@ -32,11 +32,14 @@
 //! pulls: every population of a cell comes from one source (a neighbour, a periodic image, the
 //! cell itself through a wall, or an open boundary). Cells whose 18 neighbours are all fluid
 //! cells inside the box stream by one shifted copy per direction; the others follow per-link
-//! rules. Every loop runs in a fixed order with no parallel reductions, so the result is the
-//! same bit for bit on every platform.
+//! rules. Collision and streaming run by fixed blocks of cells in parallel, each cell writing
+//! only its own populations; the wall forces are summed in link order. There are no parallel
+//! reductions, so the result is the same bit for bit on any number of threads and platform.
 
 use modelcraft_model::{BoxFace, FaceCondition, Profile};
 
+use modelcraft_parallel as parallel;
+
 use crate::lattice::{C, Q, W, c, equilibrium, moving_wall, opposite};
 use crate::setup::{MAX_LATTICE_SPEED, Obstacle, Setup};
 use crate::voxel::{Grid, TriangleTree, voxelise};
@@ -119,6 +122,24 @@ pub struct Flow {
 
 const NONE: u32 = u32::MAX;
 
+/// Cells per parallel block of the collision and the streaming. The blocks are fixed, every
+/// cell's update reads and writes only its own populations (collision) or writes only its own
+/// (streaming), so the result does not depend on the number of threads.
+const CELLS: usize = 4096;
+
+/// The first `n` populations of each direction of an array with direction stride `stride`, cut
+/// into blocks of [`CELLS`] cells: block b holds the Q directions' slices of cells
+/// b·CELLS.. .
+fn cell_blocks(data: &mut [f64], stride: usize, n: usize) -> Vec<Vec<&mut [f64]>> {
+    let mut blocks: Vec<Vec<&mut [f64]>> = (0..n.div_ceil(CELLS)).map(|_| Vec::with_capacity(Q)).collect();
+    for direction in data.chunks_mut(stride.max(1)).take(Q) {
+        for (block, cells) in blocks.iter_mut().zip(direction.get_mut(..n).unwrap_or_default().chunks_mut(CELLS)) {
+            block.push(cells);
+        }
+    }
+    blocks
+}
+
 /// A vector allocated without aborting when memory runs out.
 fn filled<T: Clone>(len: usize, value: T, what: &str) -> Result<Vec<T>> {
     let mut v = Vec::new();
@@ -423,24 +444,29 @@ impl Flow {
         Ok(())
     }
 
-    /// Two-relaxation-time collision with Guo forcing, in place. Returns the largest squared
-    /// lattice speed, or NaN when a density is not positive.
+    /// Two-relaxation-time collision with Guo forcing, in place, by blocks of [`CELLS`] cells
+    /// in parallel. Returns the largest squared lattice speed, or NaN when a density is not
+    /// positive.
     fn collide(&mut self) -> f64 {
         let n = self.n;
-        let s = self.stride;
         let (wp, wm) = (self.omega_plus, self.omega_minus);
         let (sp, sm) = (1.0 - 0.5 * wp, 1.0 - 0.5 * wm);
         let a = self.acceleration;
         let forced = a != [0.0; 3];
-        let mut max_u2: f64 = 0.0;
+        let fluid = &self.fluid;
+        let mut blocks = cell_blocks(&mut self.f, self.stride, n);
+        let mut maxima = vec![0.0f64; blocks.len()];
+        let mut work: Vec<(&mut Vec<&mut [f64]>, &mut f64)> = blocks.iter_mut().zip(maxima.iter_mut()).collect();
+        parallel::for_each_mut(&mut work, |b, (f, max_u2)| {
             let mut fq = [0.0; Q];
-        for cell in 0..n {
-            if !self.fluid[cell] {
+            let len = f.first().map_or(0, |d| d.len());
+            for i in 0..len {
+                if !fluid.get(b * CELLS + i).copied().unwrap_or(false) {
                     continue;
                 }
                 let (mut rho, mut j) = (0.0, [0.0; 3]);
                 for q in 0..Q {
-                let v = self.f[q * s + cell];
+                    let v = f.get(q).and_then(|d| d.get(i)).copied().unwrap_or(0.0);
                     fq[q] = v;
                     let cq = C[q];
                     rho += v;
@@ -449,20 +475,26 @@ impl Flow {
                     j[2] += v * f64::from(cq[2]);
                 }
                 if !(rho > 0.0 && rho.is_finite()) {
-                return f64::NAN;
+                    **max_u2 = f64::NAN;
+                    return;
                 }
                 let force = [rho * a[0], rho * a[1], rho * a[2]];
                 let u = [(j[0] + 0.5 * force[0]) / rho, (j[1] + 0.5 * force[1]) / rho, (j[2] + 0.5 * force[2]) / rho];
                 let uu = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
-            max_u2 = max_u2.max(uu);
+                **max_u2 = max_u2.max(uu);
                 let uf = u[0] * force[0] + u[1] * force[1] + u[2] * force[2];
+                let mut put = |q: usize, v: f64| {
+                    if let Some(x) = f.get_mut(q).and_then(|d| d.get_mut(i)) {
+                        *x = v;
+                    }
+                };
                 // Rest population: symmetric only.
                 let eq0 = W[0] * rho * (1.0 - 1.5 * uu);
                 let mut out0 = fq[0] - wp * (fq[0] - eq0);
                 if forced {
                     out0 += sp * W[0] * (-3.0 * uf);
                 }
-            self.f[cell] = out0;
+                put(0, out0);
                 for q in (1..Q).step_by(2) {
                     let o = q + 1;
                     let cq = C[q];
@@ -481,44 +513,60 @@ impl Flow {
                         s_plus = sp * w * (9.0 * cu * cf - 3.0 * uf);
                         s_minus = sm * w * 3.0 * cf;
                     }
-                self.f[q * s + cell] = fq[q] - relax_plus - relax_minus + s_plus + s_minus;
-                self.f[o * s + cell] = fq[o] - relax_plus + relax_minus + s_plus - s_minus;
+                    put(q, fq[q] - relax_plus - relax_minus + s_plus + s_minus);
+                    put(o, fq[o] - relax_plus + relax_minus + s_plus - s_minus);
                 }
             }
-        max_u2
+        });
+        // The largest speed of all blocks (a maximum does not depend on the order), or NaN.
+        if maxima.iter().any(|m| m.is_nan()) { f64::NAN } else { maxima.iter().fold(0.0, |m, x| m.max(*x)) }
     }
 
-    /// Pull streaming into `g`, with the wall rules and the momentum exchange.
+    /// Pull streaming into `g`, with the wall rules (blocks of [`CELLS`] destination cells in
+    /// parallel), then the momentum exchange (in link order).
     fn stream(&mut self) {
         let n = self.n;
         let s = self.stride;
+        let (f, links, offset) = (&self.f, &self.links, &self.offset);
+        let mut blocks = cell_blocks(&mut self.g, s, n);
+        parallel::for_each_mut(&mut blocks, |b, g| {
+            let start = b * CELLS;
             // Every direction as one shifted copy: right for cells whose neighbours are all fluid
             // cells in the box; the others are overwritten by their rules below.
-        for q in 0..Q {
-            let off = self.offset[q];
-            let (src, dst) = (&self.f[q * s..q * s + n], &mut self.g[q * s..q * s + n]);
-            if off >= 0 {
-                let o = off as usize;
-                if o < n {
-                    dst[o..].copy_from_slice(&src[..n - o]);
-                }
-            } else {
-                let o = off.unsigned_abs();
-                if o < n {
-                    dst[..n - o].copy_from_slice(&src[o..]);
+            for (q, dst) in g.iter_mut().enumerate() {
+                let off = offset.get(q).copied().unwrap_or(0);
+                let src = f.get(q * s..q * s + n).unwrap_or_default();
+                for (i, x) in dst.iter_mut().enumerate() {
+                    if let Some(from) = (start + i).checked_add_signed(-off)
+                        && let Some(v) = src.get(from)
+                    {
+                        *x = *v;
                     }
                 }
             }
+            let first = links.partition_point(|l| (l.cell as usize) < start);
+            for link in links.get(first..).unwrap_or_default().iter().take_while(|l| (l.cell as usize) < start + CELLS) {
+                let (cell, q) = (link.cell as usize, link.q as usize);
+                let v = match link.rule {
+                    Rule::Pull { from } => f.get(q * s + from as usize).copied(),
+                    Rule::Bounce { du, .. } => f.get(opposite(q) * s + cell).map(|out| out + du),
+                    Rule::Reflect { from, q: q2, .. } => f.get(q2 as usize * s + from as usize).copied(),
+                    // A placeholder the open boundary replaces.
+                    Rule::Open => f.get(opposite(q) * s + cell).copied(),
+                };
+                if let (Some(v), Some(x)) = (v, g.get_mut(q).and_then(|d| d.get_mut(cell - start))) {
+                    *x = v;
+                }
+            }
+        });
         self.forces.fill([0.0; 3]);
         for link in &self.links {
             let (cell, q) = (link.cell as usize, link.q as usize);
             let cq = c(q);
             match link.rule {
-                Rule::Pull { from } => self.g[q * s + cell] = self.f[q * s + from as usize],
                 Rule::Bounce { group, du } => {
                     let out = self.f[opposite(q) * s + cell];
                     let back = out + du;
-                    self.g[q * s + cell] = back;
                     if let Some(force) = self.forces.get_mut(group as usize) {
                         // c_out f* − c_in f' with c_out = −c_q.
                         for a in 0..3 {
@@ -529,7 +577,6 @@ impl Flow {
                 Rule::Reflect { from, q: q2, group } => {
                     let q2 = q2 as usize;
                     let out = self.f[q2 * s + from as usize];
-                    self.g[q * s + cell] = out;
                     let c2 = c(q2);
                     if let Some(force) = self.forces.get_mut(group as usize) {
                         for a in 0..3 {
@@ -537,8 +584,7 @@ impl Flow {
                         }
                     }
                 }
-                // A placeholder the open boundary replaces.
-                Rule::Open => self.g[q * s + cell] = self.f[opposite(q) * s + cell],
+                Rule::Pull { .. } | Rule::Open => {}
             }
         }
     }

crates/lbm/src/lib.rs

@@ -39,8 +39,10 @@
 //! the book.
 //!
 //! # Determinism
-//! Double precision, fixed loop orders, no threads and no parallel reductions: the same model
-//! gives the same result bit for bit on every platform, wasm included.
+//! Double precision and fixed loop orders. Collision and streaming run on threads by fixed blocks
+//! of cells ([`modelcraft_parallel`]), each cell's update independent of the others; the wall
+//! forces and every other sum run in a fixed order on one thread. The same model gives the same
+//! result bit for bit on any number of threads and on every platform, wasm included.
 #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)]
 #![forbid(unsafe_code)]
 // Population arrays are indexed by cell and direction counts fixed when they are allocated.

crates/parallel/Cargo.toml

@@ -0,0 +1,13 @@
+[package]
+name = "modelcraft-parallel"
+description = "ModelCraft's deterministic data parallelism: fixed work partitions on a thread pool, results independent of the thread count"
+version.workspace = true
+edition.workspace = true
+license.workspace = true
+
+[lints]
+workspace = true
+
+# The browser build runs everything on the calling thread (same results, no threads).
+[target.'cfg(not(target_arch = "wasm32"))'.dependencies]
+rayon = "1.12"

crates/parallel/src/lib.rs

@@ -0,0 +1,264 @@
+//! Deterministic data parallelism for ModelCraft's CPU solvers.
+//!
+//! The solvers promise the same result bit for bit on every run, platform and thread count
+//! (AGENTS.md, Determinism). Threads keep that promise when every floating-point value is
+//! computed by exactly one task in an order that does not depend on how many threads run them:
+//!
+//! - **Fixed partitions.** Work is split into blocks whose sizes are constants of the algorithm
+//!   (never derived from the number of threads), and each block is computed by one task.
+//! - **Ordered results.** [`map`] returns the blocks' results in block order; whatever combines
+//!   them (sums, minima, scatters into a global vector) does so on the calling thread in that
+//!   order, or in the order the sequential algorithm used.
+//! - **No shared accumulators.** Tasks write only into the disjoint slices they are handed
+//!   ([`for_each_mut`], [`chunks_mut`]); there are no atomics on floating-point values and no
+//!   reductions whose shape depends on scheduling (work stealing decides only *when* a block
+//!   runs, never *what* it computes).
+//!
+//! Under those rules a run on one thread and a run on many execute the same operations on the
+//! same operands, so the results are identical; the browser build (`wasm32`) runs the same
+//! blocks one after another on the calling thread.
+//!
+//! **Thread count.** [`with_threads`] runs a closure on a pool of that many threads (pools are
+//! cached; at most [`MAX_THREADS`]). Outside it, the helpers use the default pool:
+//! [`ENV`] (`MODELCRAFT_THREADS`) when it holds a whole number from 1 to [`MAX_THREADS`], else
+//! the available cores.
+//!
+//! Rayon (MIT OR Apache-2.0) provides the pool (docs/adr/0001-libraries.md).
+#![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)]
+#![forbid(unsafe_code)]
+
+/// Most threads a solve may use.
+pub const MAX_THREADS: usize = 256;
+
+/// The environment variable that sets the default thread count.
+pub const ENV: &str = "MODELCRAFT_THREADS";
+
+/// A thread count from text: a whole number from 1 to [`MAX_THREADS`].
+pub fn parse_threads(text: &str) -> Option<usize> {
+    text.trim().parse::<usize>().ok().filter(|n| (1..=MAX_THREADS).contains(n))
+}
+
+/// The cores this process may use (at most [`MAX_THREADS`]; 1 in the browser).
+pub fn available() -> usize {
+    #[cfg(not(target_arch = "wasm32"))]
+    {
+        std::thread::available_parallelism().map_or(1, |n| n.get()).clamp(1, MAX_THREADS)
+    }
+    #[cfg(target_arch = "wasm32")]
+    {
+        1
+    }
+}
+
+/// The default thread count: [`ENV`] when it is valid, else [`available`].
+pub fn default_threads() -> usize {
+    #[cfg(not(target_arch = "wasm32"))]
+    {
+        std::env::var(ENV).ok().and_then(|v| parse_threads(&v)).unwrap_or_else(available)
+    }
+    #[cfg(target_arch = "wasm32")]
+    {
+        1
+    }
+}
+
+/// Run `f` with `threads` worker threads (clamped to 1..=[`MAX_THREADS`]) serving the helpers of
+/// this crate. The result does not depend on `threads`; only the time does.
+pub fn with_threads<R: Send>(threads: usize, f: impl FnOnce() -> R + Send) -> R {
+    pool::with_threads(threads.clamp(1, MAX_THREADS), f)
+}
+
+/// The number of threads the helpers would use here.
+pub fn current_threads() -> usize {
+    pool::current_threads()
+}
+
+/// `f(i)` for `i` in `0..count`, run in parallel, results in index order.
+pub fn map<T: Send>(count: usize, f: impl Fn(usize) -> T + Sync + Send) -> Vec<T> {
+    pool::map(count, f)
+}
+
+/// `f(i)` for `i` in `0..count`, results in index order, run in parallel in blocks of `grain`
+/// consecutive indices (one task per block: for cheap `f`).
+pub fn map_grained<T: Send>(count: usize, grain: usize, f: impl Fn(usize) -> T + Sync + Send) -> Vec<T> {
+    let ranges: Vec<std::ops::Range<usize>> = blocks(count, grain).collect();
+    let parts = map(ranges.len(), |k| ranges.get(k).cloned().unwrap_or(0..0).map(&f).collect::<Vec<T>>());
+    parts.into_iter().flatten().collect()
+}
+
+/// `f(i, item)` for every item, run in parallel; each task gets its own item mutably (use it to
+/// hand tasks disjoint parts of the output).
+pub fn for_each_mut<T: Send>(items: &mut [T], f: impl Fn(usize, &mut T) + Sync + Send) {
+    pool::for_each_mut(items, f);
+}
+
+/// `f(k, chunk)` for the consecutive chunks of `chunk` items of `data` (the last may be
+/// shorter), run in parallel.
+pub fn chunks_mut<T: Send>(data: &mut [T], chunk: usize, f: impl Fn(usize, &mut [T]) + Sync + Send) {
+    let mut parts: Vec<&mut [T]> = data.chunks_mut(chunk.max(1)).collect();
+    for_each_mut(&mut parts, |k, part| f(k, part));
+}
+
+/// The ranges of `len` items in blocks of `block` (the last may be shorter).
+pub fn blocks(len: usize, block: usize) -> impl Iterator<Item = std::ops::Range<usize>> + Clone {
+    let block = block.max(1);
+    (0..len.div_ceil(block)).map(move |k| k * block..((k + 1) * block).min(len))
+}
+
+#[cfg(not(target_arch = "wasm32"))]
+mod pool {
+    use std::sync::{Arc, Mutex, OnceLock};
+
+    use rayon::prelude::*;
+    use rayon::{ThreadPool, ThreadPoolBuilder};
+
+    /// Pools kept for reuse, most recently used last; older ones are dropped (their threads
+    /// end) so a caller cycling through thread counts cannot pile up threads.
+    const CACHED: usize = 4;
+
+    /// Worker stack size: the element kernels keep small matrices on the stack.
+    const STACK: usize = 8 << 20;
+
+    static POOLS: Mutex<Vec<(usize, Arc<ThreadPool>)>> = Mutex::new(Vec::new());
+    static DEFAULT: OnceLock<Option<Arc<ThreadPool>>> = OnceLock::new();
+
+    fn build(threads: usize) -> Option<Arc<ThreadPool>> {
+        ThreadPoolBuilder::new().num_threads(threads).stack_size(STACK).thread_name(|i| format!("modelcraft-worker-{i}")).build().ok().map(Arc::new)
+    }
+
+    fn cached(threads: usize) -> Option<Arc<ThreadPool>> {
+        let mut pools = POOLS.lock().unwrap_or_else(std::sync::PoisonError::into_inner);
+        if let Some(at) = pools.iter().position(|(n, _)| *n == threads) {
+            let entry = pools.remove(at);
+            let pool = entry.1.clone();
+            pools.push(entry);
+            return Some(pool);
+        }
+        let pool = build(threads)?;
+        if pools.len() >= CACHED {
+            pools.remove(0);
+        }
+        pools.push((threads, pool.clone()));
+        Some(pool)
+    }
+
+    fn default_pool() -> Option<Arc<ThreadPool>> {
+        DEFAULT.get_or_init(|| build(super::default_threads())).clone()
+    }
+
+    pub fn with_threads<R: Send>(threads: usize, f: impl FnOnce() -> R + Send) -> R {
+        match cached(threads) {
+            Some(pool) => pool.install(f),
+            // No threads to be had: the helpers fall back to the default pool or to the calling
+            // thread, with the same result.
+            None => f(),
+        }
+    }
+
+    /// Run `op` inside a pool: the current one, else the default one, else on this thread.
+    fn inside<R: Send>(op: impl FnOnce() -> R + Send) -> R {
+        if rayon::current_thread_index().is_some() {
+            return op();
+        }
+        match default_pool() {
+            Some(pool) => pool.install(op),
+            None => op(),
+        }
+    }
+
+    pub fn current_threads() -> usize {
+        if rayon::current_thread_index().is_some() {
+            return rayon::current_num_threads();
+        }
+        default_pool().map_or(1, |p| p.current_num_threads())
+    }
+
+    pub fn map<T: Send>(count: usize, f: impl Fn(usize) -> T + Sync + Send) -> Vec<T> {
+        if count <= 1 {
+            return (0..count).map(f).collect();
+        }
+        inside(|| (0..count).into_par_iter().with_max_len(1).map(&f).collect())
+    }
+
+    pub fn for_each_mut<T: Send>(items: &mut [T], f: impl Fn(usize, &mut T) + Sync + Send) {
+        if items.len() <= 1 {
+            items.iter_mut().enumerate().for_each(|(i, x)| f(i, x));
+            return;
+        }
+        inside(|| items.par_iter_mut().with_max_len(1).enumerate().for_each(|(i, x)| f(i, x)));
+    }
+}
+
+#[cfg(target_arch = "wasm32")]
+mod pool {
+    pub fn with_threads<R: Send>(_threads: usize, f: impl FnOnce() -> R + Send) -> R {
+        f()
+    }
+
+    pub fn current_threads() -> usize {
+        1
+    }
+
+    pub fn map<T: Send>(count: usize, f: impl Fn(usize) -> T + Sync + Send) -> Vec<T> {
+        (0..count).map(f).collect()
+    }
+
+    pub fn for_each_mut<T: Send>(items: &mut [T], f: impl Fn(usize, &mut T) + Sync + Send) {
+        items.iter_mut().enumerate().for_each(|(i, x)| f(i, x));
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// A sum whose value depends on the order of its terms in floating point: the helpers must
+    /// give the sequential answer on any number of threads.
+    fn work() -> (Vec<f64>, f64) {
+        let terms: Vec<f64> = (0..10_000).map(|i| ((i as f64) * 0.618_033_988_7).sin() * 10f64.powi(i % 17 - 8)).collect();
+        let partial = map(blocks(terms.len(), 97).count(), |k| {
+            let r = blocks(terms.len(), 97).nth(k).unwrap_or(0..0);
+            terms[r].iter().sum::<f64>()
+        });
+        let total = partial.iter().sum::<f64>();
+        let mut scaled = terms.clone();
+        chunks_mut(&mut scaled, 61, |k, c| c.iter_mut().for_each(|x| *x *= 1.0 + k as f64 * 1e-3));
+        (scaled, total)
+    }
+
+    #[test]
+    fn results_do_not_depend_on_the_thread_count() {
+        let one = with_threads(1, work);
+        for threads in [2, 3, 7, 16] {
+            let other = with_threads(threads, work);
+            assert!(one.0.iter().zip(&other.0).all(|(a, b)| a.to_bits() == b.to_bits()));
+            assert_eq!(one.1.to_bits(), other.1.to_bits(), "{threads} threads");
+        }
+        assert_eq!(work().1.to_bits(), one.1.to_bits(), "the default pool");
+    }
+
+    #[test]
+    fn thread_counts_are_capped_and_parsed() {
+        assert_eq!(with_threads(0, current_threads), 1);
+        assert_eq!(with_threads(3, current_threads), 3);
+        assert_eq!(with_threads(usize::MAX, current_threads), MAX_THREADS);
+        assert_eq!(parse_threads(" 8 "), Some(8));
+        for bad in ["0", "-1", "257", "1e3", "", "lots", "18446744073709551616"] {
+            assert_eq!(parse_threads(bad), None, "{bad}");
+        }
+        assert!((1..=MAX_THREADS).contains(&default_threads()));
+        // Nested use runs on the pool it is in.
+        assert_eq!(with_threads(2, || map(3, |_| current_threads())), vec![2, 2, 2]);
+        assert_eq!(with_threads(3, || map_grained(10, 3, |i| i * i)), (0..10).map(|i| i * i).collect::<Vec<_>>());
+        assert_eq!(blocks(10, 4).collect::<Vec<_>>(), vec![0..4, 4..8, 8..10]);
+        assert_eq!(blocks(0, 4).count(), 0);
+        assert_eq!(blocks(3, 0).count(), 3);
+    }
+
+    #[test]
+    fn many_thread_counts_do_not_pile_up_pools() {
+        for threads in 1..=12 {
+            assert_eq!(with_threads(threads, current_threads), threads);
+        }
+    }
+}

crates/solver/Cargo.toml

@@ -11,5 +11,6 @@ workspace = true
 [dependencies]
 modelcraft-model = { workspace = true }
 modelcraft-fem = { workspace = true }
+modelcraft-parallel = { workspace = true }
 thiserror = { workspace = true }
 serde = { workspace = true }

crates/solver/src/analysis.rs

@@ -9,6 +9,8 @@ use modelcraft_fem::{DofKind, DofLayout, PointStress, von_mises};
 use modelcraft_model::{ElementKind, Model, Procedure, Resolved, Step};
 use serde::Serialize;
 
+use modelcraft_parallel as parallel;
+
 use crate::cholesky::{Factor, FactorError, UpperColumns};
 use crate::constraints::{Budget, Reduction, coupling_rows};
 use crate::eigen::{Operator, Wanted, lanczos};
@@ -164,21 +166,46 @@ struct Assembled {
     m: Option<Csr>,
 }
 
+/// Elements whose matrices are computed together before they are added.
+const ROUND: usize = 4096;
+
+/// Elements per parallel task.
+const GRAIN: usize = 16;
+
 /// Assemble the stiffness (and the mass when `mass` is set) on the mesh's node graph.
 fn assemble(mesh: &Mesh, graph: &NodeGraph, mass: Option<MassKind>) -> Result<Assembled> {
     let layout = &mesh.layout;
     let mut k = graph.zero_matrix(layout);
     let mut m = mass.map(|_| graph.zero_matrix(layout));
-    for e in &mesh.elements {
+    // Element matrices are computed in parallel, a round at a time, and added in label order
+    // (the sums are those of a sequential loop).
+    for round in parallel::blocks(mesh.elements.len(), ROUND) {
+        let elements = mesh.elements.get(round).unwrap_or_default();
+        let matrices = parallel::map_grained(elements.len(), GRAIN, |i| -> Result<_> {
+            let Some(e) = elements.get(i) else { return Err(Error::Internal("element index".into())) };
             let kernel = e.kernel()?;
-        let kinds = kernel.dofs().node;
             let at = |err| Error::element(e.label, err);
-        add_block(&mut k, layout, &e.nodes, kinds, &kernel.stiffness(&e.d).map_err(at)?.data)?;
-        if let (Some(mm), Some(kind)) = (m.as_mut(), mass) {
+            let stiffness = kernel.stiffness(&e.d).map_err(at)?.data;
+            let element_mass = match mass {
+                None => None,
+                Some(kind) => {
                     let rho = e.density.ok_or_else(|| Error::Invalid(format!("element {} has no density", e.label)))?;
+                    Some(match kind {
+                        MassKind::Consistent => kernel.mass(rho).map_err(at)?.data,
+                        MassKind::Lumped => kernel.lumped_mass(rho).map_err(at)?,
+                    })
+                }
+            };
+            Ok((kernel.dofs().node, stiffness, element_mass))
+        });
+        for (e, computed) in elements.iter().zip(matrices) {
+            let (kinds, stiffness, element_mass) = computed?;
+            add_block(&mut k, layout, &e.nodes, kinds, &stiffness)?;
+            if let (Some(mm), Some(kind), Some(data)) = (m.as_mut(), mass, element_mass) {
                 match kind {
-                MassKind::Consistent => add_scalar(mm, layout, &e.nodes, kinds, &kernel.mass(rho).map_err(at)?.data)?,
-                MassKind::Lumped => add_diagonal(mm, layout, &e.nodes, kinds, &kernel.lumped_mass(rho).map_err(at)?)?,
+                    MassKind::Consistent => add_scalar(mm, layout, &e.nodes, kinds, &data)?,
+                    MassKind::Lumped => add_diagonal(mm, layout, &e.nodes, kinds, &data)?,
+                }
             }
         }
     }
@@ -224,6 +251,8 @@ fn factor_error(e: FactorError, mesh: &Mesh, reduction: &Reduction) -> Error {
         FactorError::TooLarge { entries, limit } => {
             Error::Limit(format!("the factorization needs {entries} entries; the limit is {limit} (refine less or raise the limit)"))
         }
+        FactorError::OutOfMemory { entries } => Error::Limit(format!("not enough memory for the factorization ({entries} entries)")),
+        FactorError::Internal(what) => Error::Internal(what),
         FactorError::NotPositiveDefinite { column } => {
             let at = reduction.eq_dof.get(column).map_or_else(|| "an unknown degree of freedom".to_string(), |d| mesh.describe_dof(*d));
             Error::Singular(format!(
@@ -370,14 +399,18 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul
     let mut elements = Vec::with_capacity(mesh.elements.len());
     let mut sum = vec![[0.0; 6]; mesh.labels.len()];
     let mut count = vec![0u32; mesh.labels.len()];
-    for e in &mesh.elements {
+    let all_points = parallel::map_grained(mesh.elements.len(), GRAIN, |i| -> Result<_> {
+        let Some(e) = mesh.elements.get(i) else { return Err(Error::Internal("element index".into())) };
         let kernel = e.kernel()?;
-        let at = |err| Error::element(e.label, err);
         let ue: Vec<[f64; 3]> = e.nodes.iter().map(|n| displacement.get(*n).copied().unwrap_or_default()).collect();
         let thermal = field.as_deref().and_then(|t| e.thermal_strain(t));
-        let points = kernel.stresses_with(&e.d, &ue, thermal.as_ref()).map_err(at)?;
+        let points = kernel.stresses_with(&e.d, &ue, thermal.as_ref()).map_err(|err| Error::element(e.label, err))?;
+        Ok((points, kernel.ty(), kernel.rule()))
+    });
+    for (e, computed) in mesh.elements.iter().zip(all_points) {
+        let (points, ty, rule) = computed?;
+        let at = |err| Error::element(e.label, err);
         let values: Vec<[f64; 6]> = points.iter().map(|p| p.stress).collect();
-        let (ty, rule) = (kernel.ty(), kernel.rule());
         let cached = extrapolations.iter().position(|(t, r, _)| std::ptr::eq(*t, ty) && *r == rule);
         let index = match cached {
             Some(i) => i,

crates/solver/src/cholesky.rs

@@ -1,22 +1,58 @@
-//! Sparse Cholesky factorization A = L Lᵀ of a symmetric positive definite matrix.
+//! Sparse Cholesky factorization A = L Lᵀ of a symmetric positive definite matrix: supernodal,
+//! left-looking, multithreaded, with a fixed operation order.
 //!
-//! The up-looking method of T. A. Davis, *Direct Methods for Sparse Linear Systems*, SIAM 2006,
-//! §4.1–4.8: the elimination tree (J. W. H. Liu, *The role of elimination trees in sparse
-//! factorization*, SIAM J. Matrix Anal. Appl. 11(1):134–172, 1990) gives the pattern of each row
-//! of L as a reach in the tree, the symbolic pass counts every column before any value is
-//! computed (so the factor's size is checked against a cap before memory is committed), and each
-//! row of L is a sparse triangular solve. Operations run in a fixed order, with no threads and no
-//! data-dependent reduction order, so the factor is the same bit for bit on every platform.
+//! **Symbolic.** The elimination tree (J. W. H. Liu, *The role of elimination trees in sparse
+//! factorization*, SIAM J. Matrix Anal. Appl. 11(1):134–172, 1990) and the column counts of L
+//! come from the row subtrees (T. A. Davis, *Direct Methods for Sparse Linear Systems*, SIAM 2006,
+//! §4.1–4.4), so the factor's size is checked against a cap before any memory is committed. The
+//! columns are renumbered in a postorder of the tree (an equivalent ordering: the same fill,
+//! Liu 1990), so every subtree is a contiguous range of columns, and consecutive columns with
+//! nested structures are merged into **fundamental supernodes** (J. W. H. Liu, E. G. Ng, B. W.
+//! Peyton, *On finding supernodes for sparse matrix computations*, SIAM J. Matrix Anal. Appl.
+//! 14(1):242–252, 1993): a supernode's columns share one row structure and are stored as one
+//! dense lower trapezoid, so the factor holds exactly the entries of L.
+//!
+//! **Numeric** (left-looking supernodal Cholesky: E. Ng, B. W. Peyton, *Block sparse Cholesky
+//! algorithms on advanced uniprocessor computers*, SIAM J. Sci. Comput. 14(5):1034–1056, 1993;
+//! E. Rothberg, A. Gupta, *An efficient block-oriented approach to parallel sparse Cholesky
+//! factorization*, SIAM J. Sci. Comput. 15(6):1413–1439, 1994): a supernode starts as its
+//! columns of A, receives the update L_K L_Kᵀ of every descendant supernode K whose rows reach
+//! it (K in increasing order), and is then factored densely by panels of [`PANEL`] columns: a
+//! panel is updated with the supernode's columns before it, its diagonal block is factored column
+//! by column, and the rows below are solved against that block.
+//!
+//! **Threads and determinism.** Each entry of L is computed by one task, by a sequence of
+//! operations fixed by the matrix and the constants [`TILE`] and [`PANEL`]: sums over a
+//! supernode's columns run in column order from zero, and updates are subtracted in supernode
+//! order. Disjoint subtrees of the supernodal tree are factored concurrently (a supernode reads
+//! only its descendants), and the dense work of a large supernode is split into row tiles;
+//! neither changes any entry's sequence of operations, so the factor is the same bit for bit on
+//! any number of threads and on every platform (multiplications and additions are never fused).
 //!
 //! A pivot at or below `PIVOT_TOLERANCE` times the matrix's own diagonal entry means the matrix
 //! is singular to working precision (an unconstrained rigid-body motion, a disconnected part, a
-//! zero-stiffness direction); the factorization stops and names the column.
+//! zero-stiffness direction); the factorization stops and names the column: the first such column
+//! in the postorder, whichever thread meets it.
+
+use modelcraft_parallel as parallel;
 
 use crate::{Error, Result};
 
 /// Relative pivot size below which the matrix counts as singular.
 pub const PIVOT_TOLERANCE: f64 = 1e-10;
 
+/// Rows per tile of the dense kernels (a tile's partial sums stay in registers and L1).
+pub const TILE: usize = 64;
+
+/// Columns per panel of a supernode's dense factorization.
+pub const PANEL: usize = 64;
+
+/// Dense work (multiply–adds) from which a kernel's tiles run in parallel.
+const PARALLEL_WORK: usize = 1 << 21;
+
+/// Subtrees are handed to threads in pieces of about this fraction of the total work.
+const SUBTREE_SHARE: usize = 256;
+
 /// The upper triangle of a symmetric matrix by columns: column k holds A(i, k) for i ≤ k, rows
 /// sorted.
 #[derive(Debug, Clone)]
@@ -83,8 +119,8 @@ fn etree(a: &UpperColumns) -> Vec<usize> {
 }
 
 /// The pattern of row k of L (excluding the diagonal): the reach of column k's entries in the
-/// elimination tree. Fills `stack[top..]` in topological order and returns `top`; `mark` holds
-/// `k` for visited columns.
+/// elimination tree. Fills `stack[top..]` and returns `top`; `mark` holds `k` for visited
+/// columns.
 fn ereach(a: &UpperColumns, k: usize, parent: &[usize], stack: &mut [usize], mark: &mut [usize], path: &mut Vec<usize>) -> usize {
     let n = a.n;
     let mut top = n;
@@ -115,25 +151,6 @@ fn ereach(a: &UpperColumns, k: usize, parent: &[usize], stack: &mut [usize], mar
     top
 }
 
-/// A Cholesky factor, L by columns with the diagonal first in each column.
-#[derive(Debug, Clone)]
-pub struct Factor {
-    pub n: usize,
-    ptr: Vec<usize>,
-    idx: Vec<usize>,
-    val: Vec<f64>,
-}
-
-/// Why a factorization stopped.
-#[derive(Debug, Clone, PartialEq)]
-pub enum FactorError {
-    /// The pivot of this column is not positive (relative to the diagonal).
-    NotPositiveDefinite { column: usize },
-    /// The factor would hold more entries than allowed.
-    TooLarge { entries: usize, limit: usize },
-}
-
-impl Factor {
 /// The elimination tree and the entry count of each column of L, without computing values.
 fn symbolic(a: &UpperColumns) -> (Vec<usize>, Vec<usize>) {
     let parent = etree(a);
@@ -152,80 +169,635 @@ impl Factor {
     (parent, counts)
 }
 
+/// A postorder of a forest given by parents, children in increasing order: `post[k]` is the
+/// k-th node.
+fn postorder(parent: &[usize]) -> Vec<usize> {
+    let n = parent.len();
+    // Children as linked lists in increasing order (built from the largest child down).
+    let mut head = vec![NONE; n];
+    let mut next = vec![NONE; n];
+    for j in (0..n).rev() {
+        let Some(&p) = parent.get(j) else { continue };
+        if let (Some(h), Some(x)) = (head.get(p).copied(), next.get_mut(j)) {
+            *x = h;
+            if let Some(h) = head.get_mut(p) {
+                *h = j;
+            }
+        }
+    }
+    let mut post = Vec::with_capacity(n);
+    let mut stack = Vec::new();
+    for root in (0..n).filter(|j| parent.get(*j) == Some(&NONE)) {
+        stack.push(root);
+        while let Some(&top) = stack.last() {
+            let child = head.get(top).copied().unwrap_or(NONE);
+            if child == NONE {
+                post.push(top);
+                stack.pop();
+            } else {
+                // Detach the first remaining child and descend into it.
+                if let Some(h) = head.get_mut(top) {
+                    *h = next.get(child).copied().unwrap_or(NONE);
+                }
+                stack.push(child);
+            }
+        }
+    }
+    post
+}
+
+/// Why a factorization stopped.
+#[derive(Debug, Clone, PartialEq)]
+pub enum FactorError {
+    /// The pivot of this column is not positive (relative to the diagonal).
+    NotPositiveDefinite { column: usize },
+    /// The factor would hold more entries than allowed.
+    TooLarge { entries: usize, limit: usize },
+    /// The memory for the factor could not be had.
+    OutOfMemory { entries: usize },
+    /// A broken internal invariant (a bug), reported instead of a panic.
+    Internal(String),
+}
+
+/// One supernode's update to a later one: supernode `k`'s rows from position `p1` on reach the
+/// later one, whose columns they meet at positions `p1..p2`.
+#[derive(Debug, Clone, Copy)]
+struct Update {
+    k: usize,
+    p1: usize,
+    p2: usize,
+}
+
+/// The supernodal structure (postorder numbering).
+#[derive(Debug, Clone)]
+struct Structure {
+    /// First column of each supernode, and `n` at the end.
+    first: Vec<usize>,
+    /// Row structure of each supernode (its own columns first, then the rows below), sorted.
+    rptr: Vec<usize>,
+    rows: Vec<usize>,
+}
+
+impl Structure {
+    fn cols(&self, s: usize) -> std::ops::Range<usize> {
+        self.first.get(s).copied().unwrap_or(0)..self.first.get(s + 1).copied().unwrap_or(0)
+    }
+
+    fn rows(&self, s: usize) -> &[usize] {
+        let (a, b) = (self.rptr.get(s).copied().unwrap_or(0), self.rptr.get(s + 1).copied().unwrap_or(0));
+        self.rows.get(a..b).unwrap_or_default()
+    }
+
+    /// Entries of the supernode's lower trapezoid.
+    fn size(&self, s: usize) -> usize {
+        trapezoid(self.rows(s).len(), self.cols(s).len())
+    }
+}
+
+/// Entries of an m-row, w-column lower trapezoid (column c holds rows c..m).
+fn trapezoid(m: usize, w: usize) -> usize {
+    (0..w).map(|c| m.saturating_sub(c)).sum()
+}
+
+/// Where column c of an m-row trapezoid starts.
+fn offset(m: usize, c: usize) -> usize {
+    c * m - c * c.saturating_sub(1) / 2
+}
+
+/// A dense lower trapezoid of m rows: one supernode's values, column by column, column c
+/// holding rows c..m of the supernode's row structure.
+struct Panel<'a> {
+    m: usize,
+    data: &'a [f64],
+}
+
+impl Panel<'_> {
+    /// L(r, c) for r ≥ c.
+    fn at(&self, r: usize, c: usize) -> f64 {
+        self.data.get(offset(self.m, c) + r - c).copied().unwrap_or(0.0)
+    }
+
+    /// Rows r..r+len of column c (r ≥ c).
+    fn column(&self, c: usize, r: usize, len: usize) -> &[f64] {
+        let start = offset(self.m, c) + r - c;
+        self.data.get(start..start + len).unwrap_or_default()
+    }
+}
+
+/// S(i, c) = Σ_k L(p + i, k) L(p + c, k) over the columns k < `kw` of a trapezoid, for rows i in
+/// 0..mu and columns c in 0..wu (meaningful for i ≥ c), each sum in column order from zero.
+/// Returned by row tiles: tile t holds rows t·TILE.. as `tile[c * TILE + i]`.
+fn products(l: &Panel, kw: usize, p: usize, mu: usize, wu: usize) -> Vec<Vec<f64>> {
+    // Rows p..p+wu of the kw columns, transposed so that column c's coefficients are contiguous.
+    let mut yt = vec![0.0; wu * kw];
+    for k in 0..kw {
+        for (c, y) in l.column(k, p, wu).iter().enumerate() {
+            if let Some(slot) = yt.get_mut(c * kw + k) {
+                *slot = *y;
+            }
+        }
+    }
+    let tiles = mu.div_ceil(TILE);
+    let tile = |t: usize| -> Vec<f64> {
+        let i0 = t * TILE;
+        let h = TILE.min(mu - i0);
+        // The tile's rows of every column, padded to TILE rows.
+        let mut x = vec![0.0; kw * TILE];
+        for k in 0..kw {
+            if let Some(dst) = x.get_mut(k * TILE..k * TILE + h) {
+                dst.copy_from_slice(l.column(k, p + i0, h));
+            }
+        }
+        let cols = wu.min(i0 + h);
+        let mut out = vec![0.0; cols * TILE];
+        for c in 0..cols {
+            let mut acc = [0.0f64; TILE];
+            let ys = yt.get(c * kw..(c + 1) * kw).unwrap_or_default();
+            for (xs, b) in x.chunks_exact(TILE).zip(ys) {
+                for (a, xi) in acc.iter_mut().zip(xs) {
+                    *a += xi * b;
+                }
+            }
+            if let Some(dst) = out.get_mut(c * TILE..(c + 1) * TILE) {
+                dst.copy_from_slice(&acc);
+            }
+        }
+        out
+    };
+    if mu.saturating_mul(wu).saturating_mul(kw) >= PARALLEL_WORK { parallel::map(tiles, tile) } else { (0..tiles).map(tile).collect() }
+}
+
+/// own(rel(i), rel(c)) −= S(i, c) for i ≥ c, with S in the tiles of [`products`].
+fn subtract(own: &mut [f64], m: usize, tiles: &[Vec<f64>], mu: usize, wu: usize, rel: impl Fn(usize) -> usize) {
+    for (t, tile) in tiles.iter().enumerate() {
+        let i0 = t * TILE;
+        let h = TILE.min(mu - i0);
+        for c in 0..wu.min(i0 + h) {
+            let cc = rel(c);
+            for i in i0.max(c)..i0 + h {
+                let rr = rel(i);
+                if let (Some(x), Some(s)) = (own.get_mut(offset(m, cc) + rr - cc), tile.get(c * TILE + i - i0)) {
+                    *x -= s;
+                }
+            }
+        }
+    }
+}
+
+/// Rows c1..m of the panel c0..c1 solved against its factored diagonal block, by tiles of
+/// [`TILE`] rows (`tile[c * TILE + i]`): L(r, c) −= L(r, k) L(c, k) for the panel's k < c in
+/// order, then ÷ L(c, c).
+fn solve_below(l: &Panel, c0: usize, c1: usize) -> Vec<Vec<f64>> {
+    let (m, pw) = (l.m, c1 - c0);
+    let tiles = (m - c1).div_ceil(TILE);
+    let tile = |t: usize| -> Vec<f64> {
+        let i0 = c1 + t * TILE;
+        let h = TILE.min(m - i0);
+        let mut buf = vec![0.0; pw * TILE];
+        for c in 0..pw {
+            if let Some(dst) = buf.get_mut(c * TILE..c * TILE + h) {
+                dst.copy_from_slice(l.column(c0 + c, i0, h));
+            }
+        }
+        for c in 0..pw {
+            let (done, rest) = buf.split_at_mut(c * TILE);
+            let Some(target) = rest.get_mut(..TILE) else { continue };
+            for (k, src) in done.chunks_exact(TILE).enumerate() {
+                let lck = l.at(c0 + c, c0 + k);
+                for (y, x) in target.iter_mut().zip(src) {
+                    *y -= x * lck;
+                }
+            }
+            let d = l.at(c0 + c, c0 + c);
+            for y in target.iter_mut() {
+                *y /= d;
+            }
+        }
+        buf
+    };
+    if (m - c1).saturating_mul(pw).saturating_mul(pw) >= PARALLEL_WORK { parallel::map(tiles, tile) } else { (0..tiles).map(tile).collect() }
+}
+
+/// The matrix in postorder by lower columns: column j holds A(i, j) for i ≥ j, rows sorted.
+struct Lower {
+    ptr: Vec<usize>,
+    idx: Vec<usize>,
+    val: Vec<f64>,
+    diag: Vec<f64>,
+}
+
+impl Lower {
+    fn new(a: &UpperColumns, inverse: &[usize]) -> Lower {
+        let n = a.n;
+        let new = |i: usize| inverse.get(i).copied().unwrap_or(0);
+        let mut ptr = vec![0usize; n + 1];
+        for k in 0..n {
+            for &i in a.column(k).0 {
+                if let Some(c) = ptr.get_mut(new(i).min(new(k)) + 1) {
+                    *c += 1;
+                }
+            }
+        }
+        for j in 0..n {
+            let before = ptr.get(j).copied().unwrap_or(0);
+            if let Some(c) = ptr.get_mut(j + 1) {
+                *c += before;
+            }
+        }
+        let mut pairs = vec![(0usize, 0.0f64); ptr.last().copied().unwrap_or(0)];
+        let mut next = ptr.clone();
+        let mut diag = vec![0.0; n];
+        for k in 0..n {
+            let (rows, vals) = a.column(k);
+            for (&i, &v) in rows.iter().zip(vals) {
+                let (r, c) = (new(i).max(new(k)), new(i).min(new(k)));
+                if r == c
+                    && let Some(d) = diag.get_mut(c)
+                {
+                    *d = v;
+                }
+                if let Some(slot) = next.get_mut(c) {
+                    if let Some(p) = pairs.get_mut(*slot) {
+                        *p = (r, v);
+                    }
+                    *slot += 1;
+                }
+            }
+        }
+        for j in 0..n {
+            let (s, e) = (ptr.get(j).copied().unwrap_or(0), ptr.get(j + 1).copied().unwrap_or(0));
+            if let Some(col) = pairs.get_mut(s..e) {
+                col.sort_by_key(|p| p.0);
+            }
+        }
+        let (idx, val) = pairs.into_iter().unzip();
+        Lower { ptr, idx, val, diag }
+    }
+
+    fn column(&self, j: usize) -> (&[usize], &[f64]) {
+        let (a, b) = (self.ptr.get(j).copied().unwrap_or(0), self.ptr.get(j + 1).copied().unwrap_or(0));
+        (self.idx.get(a..b).unwrap_or_default(), self.val.get(a..b).unwrap_or_default())
+    }
+}
+
+/// What a supernode's factorization reads: the matrix, the structure and the updates each
+/// supernode receives.
+struct Context<'a> {
+    a: &'a Lower,
+    s: &'a Structure,
+    updates: &'a [Vec<Update>],
+}
+
+impl Context<'_> {
+    /// Factor supernode `j` into `own`; its descendants' values are `done[k - base]`. Fails with
+    /// the first column (postorder) whose pivot is not positive.
+    fn supernode(&self, j: usize, done: &[Vec<f64>], base: usize, own: &mut [f64]) -> std::result::Result<(), FactorError> {
+        let internal = |what: &str| FactorError::Internal(format!("supernodal factorization: {what}"));
+        let rows = self.s.rows(j);
+        let cols = self.s.cols(j);
+        let (m, w, f) = (rows.len(), cols.len(), cols.start);
+        if own.len() != trapezoid(m, w) {
+            return Err(internal("supernode size"));
+        }
+        own.fill(0.0);
+        // The columns of A.
+        for c in 0..w {
+            let (ai, av) = self.a.column(f + c);
+            for (r, v) in ai.iter().zip(av) {
+                let at = rows.binary_search(r).map_err(|_| internal("an entry of A outside the structure"))?;
+                if at < c {
+                    return Err(internal("an entry above the diagonal"));
+                }
+                if let Some(x) = own.get_mut(offset(m, c) + at - c) {
+                    *x = *v;
+                }
+            }
+        }
+        // The descendants' updates, in supernode order.
+        let mut relative = Vec::new();
+        for u in self.updates.get(j).map(Vec::as_slice).unwrap_or_default() {
+            let lk = done.get(u.k.wrapping_sub(base)).ok_or_else(|| internal("an update from outside the subtree"))?;
+            let (kr, kw) = (self.s.rows(u.k), self.s.cols(u.k).len());
+            let below = kr.get(u.p1..).unwrap_or_default();
+            let (mu, wu) = (below.len(), u.p2 - u.p1);
+            // Positions of K's rows in J's structure (both sorted).
+            relative.clear();
+            let mut at = 0;
+            for r in below {
+                while rows.get(at).is_some_and(|x| x < r) {
+                    at += 1;
+                }
+                if rows.get(at) != Some(r) {
+                    return Err(internal("an update row outside the structure"));
+                }
+                relative.push(at);
+            }
+            let tiles = products(&Panel { m: kr.len(), data: lk }, kw, u.p1, mu, wu);
+            subtract(own, m, &tiles, mu, wu, |i| relative.get(i).copied().unwrap_or(0));
+        }
+        // Dense factorization by panels.
+        for c0 in (0..w).step_by(PANEL) {
+            let c1 = (c0 + PANEL).min(w);
+            let pw = c1 - c0;
+            if c0 > 0 {
+                // The panel minus its products with the columns before it.
+                let tiles = products(&Panel { m, data: own }, c0, c0, m - c0, pw);
+                subtract(own, m, &tiles, m - c0, pw, |i| c0 + i);
+            }
+            // The diagonal block, column by column.
+            for c in c0..c1 {
+                let dst = offset(m, c);
+                for k in c0..c {
+                    let src = offset(m, k) + c - k;
+                    let lck = own.get(src).copied().unwrap_or(0.0);
+                    for r in 0..c1 - c {
+                        let x = own.get(src + r).copied().unwrap_or(0.0);
+                        if let Some(y) = own.get_mut(dst + r) {
+                            *y -= x * lck;
+                        }
+                    }
+                }
+                let d = own.get(dst).copied().unwrap_or(0.0);
+                let diag = self.a.diag.get(f + c).copied().unwrap_or(0.0);
+                if !(d > PIVOT_TOLERANCE * diag.abs()) || !d.is_finite() {
+                    return Err(FactorError::NotPositiveDefinite { column: f + c });
+                }
+                let d = d.sqrt();
+                if let Some(x) = own.get_mut(dst) {
+                    *x = d;
+                }
+                for x in own.get_mut(dst + 1..dst + c1 - c).unwrap_or_default() {
+                    *x /= d;
+                }
+            }
+            // The rows below the panel, by tiles.
+            if c1 < m {
+                let solved = solve_below(&Panel { m, data: own }, c0, c1);
+                for (t, tile) in solved.iter().enumerate() {
+                    let i0 = c1 + t * TILE;
+                    let h = TILE.min(m - i0);
+                    for c in 0..pw {
+                        let start = offset(m, c0 + c) + i0 - (c0 + c);
+                        if let (Some(dst), Some(src)) = (own.get_mut(start..start + h), tile.get(c * TILE..c * TILE + h)) {
+                            dst.copy_from_slice(src);
+                        }
+                    }
+                }
+            }
+        }
+        Ok(())
+    }
+}
+
+/// A Cholesky factor: L in supernodes, for the matrix with its columns in postorder.
+#[derive(Debug, Clone)]
+pub struct Factor {
+    pub n: usize,
+    /// `perm[k]`: the original column that is column k of the factor.
+    perm: Vec<usize>,
+    structure: Structure,
+    /// The values of each supernode (a lower trapezoid by columns).
+    values: Vec<Vec<f64>>,
+    entries: usize,
+}
+
+impl Factor {
     pub fn new(a: &UpperColumns, max_entries: usize) -> std::result::Result<Factor, FactorError> {
         let n = a.n;
-        let (parent, counts) = Factor::symbolic(a);
-        let total: usize = counts.iter().sum();
+        let internal = |what: &str| FactorError::Internal(format!("supernodal structure: {what}"));
+        let (parent0, counts0) = symbolic(a);
+        let total: usize = counts0.iter().sum();
         if total > max_entries {
             return Err(FactorError::TooLarge { entries: total, limit: max_entries });
         }
-        let mut ptr = Vec::with_capacity(n + 1);
-        ptr.push(0);
-        for c in &counts {
-            ptr.push(ptr.last().copied().unwrap_or(0) + c);
+        // Postorder: the same tree and counts, renumbered.
+        let perm = postorder(&parent0);
+        if perm.len() != n {
+            return Err(internal("postorder"));
         }
-        let mut idx = vec![0usize; total];
-        let mut val = vec![0.0f64; total];
-        // Next free slot in each column.
-        let mut next: Vec<usize> = ptr.get(..n).unwrap_or_default().to_vec();
-        let mut x = vec![0.0f64; n];
-        let mut stack = vec![0usize; n];
+        let mut inverse = vec![0usize; n];
+        for (k, &j) in perm.iter().enumerate() {
+            if let Some(x) = inverse.get_mut(j) {
+                *x = k;
+            }
+        }
+        let parent: Vec<usize> =
+            perm.iter().map(|&j| parent0.get(j).copied().filter(|p| *p != NONE).and_then(|p| inverse.get(p).copied()).unwrap_or(NONE)).collect();
+        let counts: Vec<usize> = perm.iter().map(|&j| counts0.get(j).copied().unwrap_or(1)).collect();
+        let a = Lower::new(a, &inverse);
+
+        // Fundamental supernodes: column j joins column j − 1's supernode when it is j − 1's
+        // parent, j − 1 is its only child, and their structures nest.
+        let mut children = vec![0usize; n];
+        for &p in &parent {
+            if let Some(c) = children.get_mut(p) {
+                *c += 1;
+            }
+        }
+        let mut first = vec![0usize];
+        for j in 1..n {
+            let joins = parent.get(j - 1) == Some(&j) && children.get(j) == Some(&1) && counts.get(j - 1).copied() == counts.get(j).map(|c| c + 1);
+            if !joins {
+                first.push(j);
+            }
+        }
+        if n > 0 {
+            first.push(n);
+        }
+        let ns = first.len().saturating_sub(1);
+        let mut snode = vec![0usize; n];
+        for s in 0..ns {
+            for j in first.get(s).copied().unwrap_or(0)..first.get(s + 1).copied().unwrap_or(0) {
+                if let Some(x) = snode.get_mut(j) {
+                    *x = s;
+                }
+            }
+        }
+        // The supernodal tree and the row structures: a supernode's rows are its columns, the
+        // rows of A below them and its children's rows past their own columns.
+        let sparent: Vec<usize> = (0..ns)
+            .map(|s| {
+                let last = first.get(s + 1).copied().unwrap_or(1).saturating_sub(1);
+                parent.get(last).copied().filter(|p| *p != NONE).and_then(|p| snode.get(p).copied()).unwrap_or(NONE)
+            })
+            .collect();
+        let mut kids: Vec<Vec<usize>> = vec![Vec::new(); ns];
+        for (s, &p) in sparent.iter().enumerate() {
+            if let Some(k) = kids.get_mut(p) {
+                k.push(s);
+            }
+        }
+        let mut rptr = vec![0usize];
+        let mut rows: Vec<usize> = Vec::new();
         let mut mark = vec![NONE; n];
-        let mut path = Vec::new();
-        for k in 0..n {
-            let top = ereach(a, k, &parent, &mut stack, &mut mark, &mut path);
-            let (rows, vals) = a.column(k);
-            let mut diag = 0.0;
-            for (i, v) in rows.iter().zip(vals) {
-                if *i < k {
-                    if let Some(xi) = x.get_mut(*i) {
-                        *xi = *v;
+        let mut set = Vec::new();
+        for s in 0..ns {
+            let (f, l) = (first.get(s).copied().unwrap_or(0), first.get(s + 1).copied().unwrap_or(0));
+            set.clear();
+            let mut add = |r: usize, set: &mut Vec<usize>| {
+                if let Some(x) = mark.get_mut(r)
+                    && *x != s
+                {
+                    *x = s;
+                    set.push(r);
                 }
-                } else if *i == k {
-                    diag = *v;
+            };
+            for j in f..l {
+                add(j, &mut set);
+                for &r in a.column(j).0 {
+                    add(r, &mut set);
                 }
             }
-            let mut d = diag;
-            for &i in stack.get(top..).unwrap_or_default() {
-                let start = ptr.get(i).copied().unwrap_or(0);
-                let lii = val.get(start).copied().unwrap_or(1.0);
-                let lki = x.get(i).copied().unwrap_or(0.0) / lii;
-                if let Some(xi) = x.get_mut(i) {
-                    *xi = 0.0;
-                }
-                let end = next.get(i).copied().unwrap_or(start);
-                for p in start + 1..end {
-                    let (r, l) = (idx.get(p).copied().unwrap_or(0), val.get(p).copied().unwrap_or(0.0));
-                    if let Some(xr) = x.get_mut(r) {
-                        *xr -= l * lki;
+            for &c in kids.get(s).map(Vec::as_slice).unwrap_or_default() {
+                let (cs, ce) = (rptr.get(c).copied().unwrap_or(0), rptr.get(c + 1).copied().unwrap_or(0));
+                let cw = first.get(c + 1).copied().unwrap_or(0) - first.get(c).copied().unwrap_or(0);
+                for &r in rows.get(cs + cw..ce).unwrap_or_default() {
+                    add(r, &mut set);
                 }
             }
-                d -= lki * lki;
-                if let Some(slot) = next.get_mut(i) {
-                    if let (Some(ri), Some(vi)) = (idx.get_mut(*slot), val.get_mut(*slot)) {
-                        *ri = k;
-                        *vi = lki;
+            set.sort_unstable();
+            if set.len() != counts.get(f).copied().unwrap_or(0) || set.first() != Some(&f) {
+                return Err(internal("a supernode's rows differ from its column count"));
             }
-                    *slot += 1;
+            rows.extend_from_slice(&set);
+            rptr.push(rows.len());
+        }
+        let structure = Structure { first, rptr, rows };
+
+        // The updates each supernode receives, from its descendants in increasing order.
+        let mut updates: Vec<Vec<Update>> = vec![Vec::new(); ns];
+        for k in 0..ns {
+            let kr = structure.rows(k);
+            let mut p = structure.cols(k).len();
+            while let Some(target) = kr.get(p).and_then(|r| snode.get(*r)).copied() {
+                let end = structure.cols(target).end;
+                let mut q = p;
+                while kr.get(q).is_some_and(|r| *r < end) {
+                    q += 1;
+                }
+                if let Some(list) = updates.get_mut(target) {
+                    list.push(Update { k, p1: p, p2: q });
+                }
+                p = q;
             }
         }
-            if !(d > PIVOT_TOLERANCE * diag.abs()) || !d.is_finite() {
-                return Err(FactorError::NotPositiveDefinite { column: k });
+
+        // Memory for the values, taken without aborting.
+        let mut values: Vec<Vec<f64>> = Vec::new();
+        values.try_reserve_exact(ns).map_err(|_| FactorError::OutOfMemory { entries: total })?;
+        for s in 0..ns {
+            let mut v = Vec::new();
+            v.try_reserve_exact(structure.size(s)).map_err(|_| FactorError::OutOfMemory { entries: total })?;
+            v.resize(structure.size(s), 0.0);
+            values.push(v);
         }
-            if let Some(slot) = next.get_mut(k) {
-                if let (Some(ri), Some(vi)) = (idx.get_mut(*slot), val.get_mut(*slot)) {
-                    *ri = k;
-                    *vi = d.sqrt();
+
+        // Tasks: subtrees (contiguous ranges of supernodes in postorder) of at most a share of the
+        // work each; the supernodes above them are factored afterwards, in order.
+        let mut subtree: Vec<usize> = (0..ns)
+            .map(|s| {
+                let (m, w) = (structure.rows(s).len(), structure.cols(s).len());
+                (0..w).map(|c| (m - c).saturating_mul(m - c)).fold(0usize, usize::saturating_add)
+            })
+            .collect();
+        let mut start: Vec<usize> = (0..ns).collect();
+        for s in 0..ns {
+            if let Some(&p) = sparent.get(s)
+                && p != NONE
+            {
+                let (ws, ss) = (subtree.get(s).copied().unwrap_or(0), start.get(s).copied().unwrap_or(s));
+                if let Some(x) = subtree.get_mut(p) {
+                    *x = x.saturating_add(ws);
                 }
-                *slot += 1;
+                if let Some(x) = start.get_mut(p) {
+                    *x = (*x).min(ss);
                 }
             }
-        Ok(Factor { n, ptr, idx, val })
+        }
+        let all = (0..ns).filter(|s| sparent.get(*s) == Some(&NONE)).filter_map(|s| subtree.get(s).copied()).fold(0, usize::saturating_add);
+        let share = (all / SUBTREE_SHARE).max(PARALLEL_WORK);
+        let fits = |s: usize| subtree.get(s).is_some_and(|w| *w <= share);
+        let mut tasks: Vec<std::ops::Range<usize>> = Vec::new();
+        let mut top: Vec<usize> = Vec::new();
+        let mut s = 0;
+        while s < ns {
+            // The largest fitting subtree whose first supernode is s.
+            let mut root = NONE;
+            let mut r = s;
+            while r != NONE && start.get(r) == Some(&s) && fits(r) {
+                root = r;
+                r = sparent.get(r).copied().unwrap_or(NONE);
+            }
+            if root == NONE {
+                top.push(s);
+                s += 1;
+            } else {
+                tasks.push(s..root + 1);
+                s = root + 1;
+            }
+        }
+
+        let context = Context { a: &a, s: &structure, updates: &updates };
+        // The tasks run in parallel, each through its supernodes in order.
+        let mut pieces: Vec<(std::ops::Range<usize>, &mut [Vec<f64>], Option<usize>)> = Vec::with_capacity(tasks.len());
+        let mut rest: &mut [Vec<f64>] = &mut values;
+        let mut at = 0;
+        for range in &tasks {
+            let (_, tail) = std::mem::take(&mut rest).split_at_mut(range.start - at);
+            let (piece, tail) = tail.split_at_mut(range.len());
+            pieces.push((range.clone(), piece, None));
+            rest = tail;
+            at = range.end;
+        }
+        let broken = std::sync::Mutex::new(None::<FactorError>);
+        parallel::for_each_mut(&mut pieces, |_, (range, piece, singular)| {
+            for j in range.clone() {
+                let (done, own) = piece.split_at_mut(j - range.start);
+                let Some(own) = own.first_mut() else { break };
+                match context.supernode(j, done, range.start, own) {
+                    Ok(()) => {}
+                    Err(FactorError::NotPositiveDefinite { column }) => {
+                        *singular = Some(column);
+                        break;
+                    }
+                    Err(e) => {
+                        if let Ok(mut b) = broken.lock() {
+                            b.get_or_insert(e);
+                        }
+                        break;
+                    }
+                }
+            }
+        });
+        if let Some(e) = broken.into_inner().unwrap_or_else(std::sync::PoisonError::into_inner) {
+            return Err(e);
+        }
+        // The first singular column in postorder is the one a pass in order would meet first:
+        // the supernodes above the tasks run up to it.
+        let earliest = pieces.iter().filter_map(|p| p.2).min();
+        drop(pieces);
+        let original = |column: usize| perm.get(column).copied().unwrap_or(column);
+        for &j in &top {
+            if earliest.is_some_and(|c| structure.cols(j).start > c) {
+                break;
+            }
+            let (done, own) = values.split_at_mut(j);
+            let Some(own) = own.first_mut() else { break };
+            match context.supernode(j, done, 0, own) {
+                Ok(()) => {}
+                Err(FactorError::NotPositiveDefinite { column }) => return Err(FactorError::NotPositiveDefinite { column: original(column) }),
+                Err(e) => return Err(e),
+            }
+        }
+        if let Some(column) = earliest {
+            return Err(FactorError::NotPositiveDefinite { column: original(column) });
+        }
+        Ok(Factor { n, perm, structure, values, entries: total })
     }
 
     /// Entries in L.
     pub fn entries(&self) -> usize {
-        self.val.len()
+        self.entries
     }
 
     /// Solve A x = b in place.
@@ -233,30 +805,43 @@ impl Factor {
         if b.len() != self.n {
             return Err(Error::Internal("right-hand side length".into()));
         }
+        let mut x: Vec<f64> = self.perm.iter().map(|&j| b.get(j).copied().unwrap_or(0.0)).collect();
+        let s = &self.structure;
         // L y = b, by columns.
-        for j in 0..self.n {
-            let (start, end) = (self.ptr.get(j).copied().unwrap_or(0), self.ptr.get(j + 1).copied().unwrap_or(0));
-            let yj = b.get(j).copied().unwrap_or(0.0) / self.val.get(start).copied().unwrap_or(1.0);
-            if let Some(x) = b.get_mut(j) {
-                *x = yj;
+        for (k, l) in self.values.iter().enumerate() {
+            let rows = s.rows(k);
+            let m = rows.len();
+            for (c, j) in s.cols(k).enumerate() {
+                let start = offset(m, c);
+                let yj = x.get(j).copied().unwrap_or(0.0) / l.get(start).copied().unwrap_or(1.0);
+                if let Some(v) = x.get_mut(j) {
+                    *v = yj;
+                }
+                for (r, v) in rows.get(c + 1..).unwrap_or_default().iter().zip(l.get(start + 1..).unwrap_or_default()) {
+                    if let Some(t) = x.get_mut(*r) {
+                        *t -= v * yj;
                     }
-            for p in start + 1..end {
-                let (r, l) = (self.idx.get(p).copied().unwrap_or(0), self.val.get(p).copied().unwrap_or(0.0));
-                if let Some(x) = b.get_mut(r) {
-                    *x -= l * yj;
                 }
             }
         }
         // Lᵀ x = y.
-        for j in (0..self.n).rev() {
-            let (start, end) = (self.ptr.get(j).copied().unwrap_or(0), self.ptr.get(j + 1).copied().unwrap_or(0));
-            let mut s = b.get(j).copied().unwrap_or(0.0);
-            for p in start + 1..end {
-                let (r, l) = (self.idx.get(p).copied().unwrap_or(0), self.val.get(p).copied().unwrap_or(0.0));
-                s -= l * b.get(r).copied().unwrap_or(0.0);
+        for (k, l) in self.values.iter().enumerate().rev() {
+            let rows = s.rows(k);
+            let m = rows.len();
+            for (c, j) in s.cols(k).enumerate().rev() {
+                let start = offset(m, c);
+                let mut sum = x.get(j).copied().unwrap_or(0.0);
+                for (r, v) in rows.get(c + 1..).unwrap_or_default().iter().zip(l.get(start + 1..).unwrap_or_default()) {
+                    sum -= v * x.get(*r).copied().unwrap_or(0.0);
                 }
-            if let Some(x) = b.get_mut(j) {
-                *x = s / self.val.get(start).copied().unwrap_or(1.0);
+                if let Some(t) = x.get_mut(j) {
+                    *t = sum / l.get(start).copied().unwrap_or(1.0);
+                }
+            }
+        }
+        for (&j, v) in self.perm.iter().zip(&x) {
+            if let Some(t) = b.get_mut(j) {
+                *t = *v;
             }
         }
         Ok(())
@@ -285,6 +870,16 @@ mod tests {
         Csr { n, ptr, idx, val }
     }
 
+    fn check_solve(csr: &Csr, f: &Factor, tol: f64) {
+        let n = csr.n;
+        let x: Vec<f64> = (0..n).map(|i| (i as f64 * 0.37).sin()).collect();
+        let mut b = vec![0.0; n];
+        csr.mul(&x, &mut b);
+        f.solve(&mut b).unwrap();
+        let err = b.iter().zip(&x).map(|(p, q)| (p - q).abs()).fold(0.0, f64::max);
+        assert!(err < tol, "{err}");
+    }
+
     #[test]
     fn factor_and_solve_a_sparse_spd_matrix() {
         // A 2D Laplacian-like matrix on a 7 × 5 grid plus a diagonal shift: sparse with fill.
@@ -307,11 +902,7 @@ mod tests {
         }
         let csr = dense_to_csr(&a);
         let f = Factor::new(&UpperColumns::from_symmetric(&csr), usize::MAX).unwrap();
-        let x: Vec<f64> = (0..n).map(|i| (i as f64 * 0.37).sin()).collect();
-        let mut b = vec![0.0; n];
-        csr.mul(&x, &mut b);
-        f.solve(&mut b).unwrap();
-        assert!(b.iter().zip(&x).all(|(p, q)| (p - q).abs() < 1e-13), "{b:?}");
+        check_solve(&csr, &f, 1e-13);
         // The cap is checked before the numeric factorization.
         assert!(matches!(Factor::new(&UpperColumns::from_symmetric(&csr), 10), Err(FactorError::TooLarge { .. })));
         // A singular matrix (a pure Laplacian has the constant vector in its null space).
@@ -322,5 +913,104 @@ mod tests {
         }
         let e = Factor::new(&UpperColumns::from_symmetric(&dense_to_csr(&sing)), usize::MAX).unwrap_err();
         assert_eq!(e, FactorError::NotPositiveDefinite { column: n - 1 });
+        // Empty and 1 × 1 matrices.
+        let empty = Csr { n: 0, ptr: vec![0], idx: vec![], val: vec![] };
+        let f = Factor::new(&UpperColumns::from_symmetric(&empty), 10).unwrap();
+        assert_eq!(f.entries(), 0);
+        f.solve(&mut []).unwrap();
+        let f = Factor::new(&UpperColumns::from_symmetric(&dense_to_csr(&[vec![4.0]])), 10).unwrap();
+        let mut b = [2.0];
+        f.solve(&mut b).unwrap();
+        assert_eq!(b, [0.5]);
+    }
+
+    /// Two 3D grids (3 coupled unknowns per point, in grid order) joined by a dense separator
+    /// ordered last, as nested dissection orders them, and a separate 5 × 5 block: subtrees run
+    /// as parallel tasks, the separator is a supernode of several panels whose kernels run in
+    /// parallel tiles. The factor solves the system and is the same bit for bit on any number of
+    /// threads.
+    #[test]
+    fn supernodes_panels_and_threads() {
+        let (nx, ny, nz) = (8, 8, 6);
+        let point = |i: usize, j: usize, k: usize| i + nx * (j + ny * k);
+        let per = 3 * nx * ny * nz;
+        let sep = 400;
+        let n = 2 * per + sep + 5;
+        let mut rows: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];
+        let pair = |rows: &mut Vec<Vec<(usize, f64)>>, a: usize, b: usize, v: f64| {
+            rows[a].push((b, v));
+            rows[b].push((a, v));
+        };
+        for base in [0, per] {
+            for k in 0..nz {
+                for j in 0..ny {
+                    for i in 0..nx {
+                        let p = base + 3 * point(i, j, k);
+                        for a in 0..3 {
+                            rows[p + a].push((p + a, 7.0));
+                            for b in a + 1..3 {
+                                pair(&mut rows, p + a, p + b, 0.25);
+                            }
+                        }
+                        let next = [
+                            (i + 1 < nx).then(|| point(i + 1, j, k)),
+                            (j + 1 < ny).then(|| point(i, j + 1, k)),
+                            (k + 1 < nz).then(|| point(i, j, k + 1)),
+                        ];
+                        for q in next.into_iter().flatten() {
+                            for a in 0..3 {
+                                pair(&mut rows, p + a, base + 3 * q + a, -1.0);
+                            }
+                        }
+                    }
+                }
+            }
+        }
+        let s0 = 2 * per;
+        for a in 0..sep {
+            rows[s0 + a].push((s0 + a, 20.0));
+            for b in a + 1..sep {
+                pair(&mut rows, s0 + a, s0 + b, 0.5 / (1.0 + (b - a) as f64).powi(2));
+            }
+            pair(&mut rows, s0 + a, (7 * a) % per, -0.1);
+            pair(&mut rows, s0 + a, per + (11 * a) % per, -0.1);
+        }
+        let d0 = s0 + sep;
+        for a in 0..5 {
+            rows[d0 + a].push((d0 + a, 10.0));
+            for b in a + 1..5 {
+                pair(&mut rows, d0 + a, d0 + b, 1.0 / (1.0 + (a + b) as f64));
+            }
+        }
+        let mut ptr = vec![0];
+        let mut idx = Vec::new();
+        let mut val = Vec::new();
+        for row in &mut rows {
+            row.sort_by_key(|e| e.0);
+            for (c, v) in row.iter() {
+                idx.push(*c);
+                val.push(*v);
+            }
+            ptr.push(idx.len());
+        }
+        let csr = Csr { n, ptr, idx, val };
+        let upper = UpperColumns::from_symmetric(&csr);
+        let one = parallel::with_threads(1, || Factor::new(&upper, usize::MAX)).unwrap();
+        let supernodes = one.structure.first.len() - 1;
+        let widest = (0..supernodes).map(|s| one.structure.cols(s).len()).max().unwrap();
+        assert!(widest > 2 * PANEL, "a supernode spans panels: {widest}");
+        check_solve(&csr, &one, 1e-12);
+        for threads in [2, 7] {
+            let other = parallel::with_threads(threads, || Factor::new(&upper, usize::MAX)).unwrap();
+            assert!(one.values.iter().flatten().zip(other.values.iter().flatten()).all(|(a, b)| a.to_bits() == b.to_bits()), "{threads} threads");
+        }
+    }
+
+    #[test]
+    fn postorder_lists_children_first() {
+        assert_eq!(postorder(&[2, 2, 4, 4, NONE, NONE]), vec![0, 1, 2, 3, 4, 5]);
+        assert_eq!(postorder(&[4, 3, 3, 4, NONE]), vec![0, 1, 2, 3, 4]);
+        assert_eq!(postorder(&[NONE, 0]), vec![1, 0]);
+        assert_eq!(postorder(&[3, NONE, 1, NONE]), vec![2, 1, 0, 3]);
     }
 }

crates/solver/src/eigen.rs

@@ -184,7 +184,7 @@ pub struct Wanted {
 }
 
 /// The operators the Lanczos iteration needs.
-pub trait Operator {
+pub trait Operator: Sync {
     fn n(&self) -> usize;
     /// y = M x.
     fn mass(&self, x: &[f64], y: &mut [f64]);
@@ -285,10 +285,12 @@ pub fn lanczos(op: &impl Operator, shift: f64, wanted: Wanted, block: usize, max
         let j0 = q.len();
         // W = Op Q_j.
         let mut w: Vec<Vec<f64>> = Vec::with_capacity(b);
-        for mx in &mcur {
-            let mut y = mx.clone();
-            op.solve(&mut y)?;
-            w.push(y);
+        // The block's solves are independent: one task each.
+        for y in modelcraft_parallel::map(mcur.len(), |i| {
+            let mut y = mcur.get(i).cloned().unwrap_or_default();
+            op.solve(&mut y).map(|()| y)
+        }) {
+            w.push(y?);
         }
         // A_j = Q_jᵀ M W (symmetrized), W −= Q_j A_j (the earlier blocks go in the full
         // reorthogonalization below).

crates/solver/src/heat.rs

@@ -308,6 +308,8 @@ fn factor(a: &Csr, limits: &Limits, mesh: &HeatMesh, eqs: &Equations) -> Result<
         FactorError::TooLarge { entries, limit } => {
             Error::Limit(format!("the factorization needs {entries} entries; the limit is {limit} (refine less or raise the limit)"))
         }
+        FactorError::OutOfMemory { entries } => Error::Limit(format!("not enough memory for the factorization ({entries} entries)")),
+        FactorError::Internal(what) => Error::Internal(what),
         FactorError::NotPositiveDefinite { column } => {
             let label = eqs.node.get(column).and_then(|n| mesh.labels.get(*n)).copied().unwrap_or(0);
             Error::Singular(format!(

crates/solver/src/lib.rs

@@ -15,8 +15,11 @@
 //!    von Mises. Frequency analysis ([`eigen`]) runs shift-invert block Lanczos on the same
 //!    kind of factorization.
 //!
-//! Everything is sequential and in a fixed order, so the same model gives the same numbers bit
-//! for bit; sizes and iterations are capped by [`Limits`].
+//! Element matrices and stresses, the factorization's subtrees and dense tiles, the Lanczos
+//! block solves and matrix products run on threads ([`modelcraft_parallel`]), but every number
+//! is computed by one task in an order fixed by the model alone and sums are taken in a fixed
+//! order, so the same model gives the same numbers bit for bit on any number of threads and on
+//! every platform; sizes and iterations are capped by [`Limits`].
 #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)]
 #![forbid(unsafe_code)]
 // Matrix kernels read as the textbook formulas with index loops.

crates/solver/src/sparse.rs

@@ -11,6 +11,9 @@ use modelcraft_fem::{DofLayout, DofSet};
 
 use crate::{Error, Result};
 
+/// Rows per parallel block of a product.
+const ROWS: usize = 4096;
+
 /// A square sparse matrix in compressed rows with sorted column indices.
 #[derive(Debug, Clone, PartialEq)]
 pub struct Csr {
@@ -28,11 +31,15 @@ impl Csr {
     }
 
     /// y = A x.
+    /// Rows are computed in parallel blocks; each row's sum runs in column order.
     pub fn mul(&self, x: &[f64], y: &mut [f64]) {
-        for (r, yr) in y.iter_mut().enumerate().take(self.n) {
-            let (cols, vals) = self.row(r);
+        let n = self.n.min(y.len());
+        modelcraft_parallel::chunks_mut(y.get_mut(..n).unwrap_or_default(), ROWS, |b, block| {
+            for (i, yr) in block.iter_mut().enumerate() {
+                let (cols, vals) = self.row(b * ROWS + i);
                 *yr = cols.iter().zip(vals).map(|(c, v)| v * x.get(*c).copied().unwrap_or(0.0)).sum();
             }
+        });
     }
 
     pub fn diagonal(&self) -> Vec<f64> {

docs/adr/0001-libraries.md

@@ -38,7 +38,28 @@ build. Our code is sequential with a fixed operation order, `unsafe`-free and wa
 factor size is known (and capped) before any value is computed. Cost if this is wrong: large
 models (beyond a few hundred thousand equations) factor slower than they could; the remedy is a
 supernodal factorization of our own, or `faer` behind the same `Factor` interface with a stated
-tolerance instead of bit-identity, as the GPU exception already does.
+tolerance instead of bit-identity, as the GPU exception already does. *Updated 2026-10-13 (CPU
+parallelism task):* the factorization is now our own supernodal left-looking Cholesky (Liu–Ng–Peyton
+1993 fundamental supernodes on a postordered elimination tree, Ng–Peyton 1993, Rothberg–Gupta
+1994), multithreaded with a fixed operation order (below). Its rounding differs from the old
+up-looking factor in the last digits; every oracle fixture stays within its band.
+
+**Threads: rayon** (MIT OR Apache-2.0), in `modelcraft-parallel` only (an L0 crate with no
+workspace dependencies; *decided 2026-10-13*). The solvers keep the bit-for-bit rule on any
+number of threads by construction: work is cut into blocks whose sizes are constants of the
+algorithm, never derived from the thread count; each floating-point value is computed by one task
+in an order fixed by the input; results come back in block order and every sum over blocks or
+elements (nodal forces, assembly, wall forces, energies) runs on one thread in the order the
+sequential algorithm used. No float atomics, no `reduce` whose tree depends on work stealing:
+rayon decides only when a block runs. `std::thread::scope` was the alternative and needs no
+dependency, but it spawns threads per parallel region; the explicit solver enters several regions
+per time increment (10⁵–10⁷ increments), so a persistent pool is needed, and writing our own would
+re-implement rayon's. rayon was already in the shipped dependency tree (THIRD-PARTY-LICENSES.txt
+unchanged). In the browser build (`wasm32`) the crate does not depend on rayon and runs the same
+blocks one after another, with the same results. The thread count is a session and `solve.*`
+parameter (`threads`, 1–256) whose default is `MODELCRAFT_THREADS` or the available cores. Cost if
+this is wrong: rayon's pool is a process-wide resource and nested pools add latency; replacing it
+touches only `modelcraft-parallel`.
 
 **Allowed:** MIT, BSD-2/3, Apache-2.0, ISC, Zlib, Unicode, MPL-2.0 (as an unmodified dependency),
 CC0 and public domain.

docs/architecture.md

@@ -19,6 +19,7 @@ crates in strictly lower layers.
 | Layer | Crate | Role | Status |
 |---|---|---|---|
 | L0 | `modelcraft-model` | Solver-neutral FE model: nodes and elements by label, their geometry tags, sets, materials, sections, couplings, BCs, loads and steps, the model check and the resolved view the writers read. Knows no file format and no geometry type. | started |
+| L0 | `modelcraft-parallel` | Deterministic data parallelism for the solvers: fixed work partitions on a rayon pool, results in block order, the thread-count setting (`MODELCRAFT_THREADS`, `threads`); sequential in the browser. No workspace dependencies. | started |
 | L1 | `modelcraft-geometry` | The only door to SolidCraft: opens `.solidcraft` designs (rebuilt from their timeline, bodies with persistent face names), B-rep bodies, measures, tessellation into a face-tagged `SurfaceMesh`; face-tagged 3MF files read into the same type without the kernel (docs/formats/tagged-mesh.md). No SolidCraft type in its API. | started |
 | L1 | `modelcraft-mesh` | Surface preparation (weld, watertightness, feature edges), robust predicates, Delaunay tetrahedralization and refinement, Voronoi-based optimization, hex meshing, quality metrics. | started |
 | L1 | `modelcraft-fem` | Element kernels shared by every structural solver: the element table (tet4, tet10, wedge6, hex8 full and incompatible modes, hex20 full and reduced), shape functions, quadrature, element stiffness, consistent and lumped mass, body and pressure loads, stresses and their extrapolation, isotropic and orthotropic elasticity, the uniform-strain (one-point) hex8 and tet4 kinematics, hourglass shape vectors and element frequency bounds of the explicit solver, objective stress rates and J2/Johnson–Cook return mapping. Pure functions, no solver loop. | started |
@@ -185,8 +186,9 @@ module docs; ADR 0001 records why the linear algebra is our own.
    symmetric positive definite (no penalty, no multipliers). Free mesh equations follow the
    node order, reference point equations come last. Without couplings T is a selection and the
    numbers are those of plain elimination, bit for bit.
-4. **Static solution.** Up-looking sparse Cholesky (the factor size is computed and capped before
-   any value), reactions T_pᵀ (K u − f) (at nodes and at prescribed reference points, forces and
+4. **Static solution.** Supernodal left-looking sparse Cholesky on the postordered elimination
+   tree (the factor size is computed and capped before any value; subtrees and dense row tiles on
+   threads with a fixed operation order, so the factor is the same on any thread count), reactions T_pᵀ (K u − f) (at nodes and at prescribed reference points, forces and
    moments), integration-point stresses extrapolated to the nodes and averaged, von Mises. Shells
    report section forces and bottom, middle and top stresses at their nodes (their top surface
    joins the nodal average), beams their end section forces; nodes report rotations and reaction

docs/solvers/literature.md

@@ -974,7 +974,10 @@ the module docs):
 - Ordering: A. George, SIAM J. Numer. Anal. 10(2):345–363, 1973; H. D. Simon, Computing Systems
   in Engineering 2:135–148, 1991.
 - Sparse Cholesky: T. A. Davis, *Direct Methods for Sparse Linear Systems*, SIAM 2006; J. W. H.
-  Liu, SIAM J. Matrix Anal. Appl. 11(1):134–172, 1990.
+  Liu, SIAM J. Matrix Anal. Appl. 11(1):134–172, 1990; supernodes: J. W. H. Liu, E. G. Ng, B. W.
+  Peyton, SIAM J. Matrix Anal. Appl. 14(1):242–252, 1993; left-looking supernodal factorization:
+  E. Ng, B. W. Peyton, SIAM J. Sci. Comput. 14(5):1034–1056, 1993; parallel block factorization:
+  E. Rothberg, A. Gupta, SIAM J. Sci. Comput. 15(6):1413–1439, 1994.
 - Eigenvalues: T. Ericsson, A. Ruhe, Math. Comp. 35(152):1251–1268, 1980; R. G. Grimes, J. G.
   Lewis, H. D. Simon, SIAM J. Matrix Anal. Appl. 15(1):228–272, 1994; B. N. Parlett, *The
   Symmetric Eigenvalue Problem*, SIAM 1998; G. H. Golub, C. F. Van Loan, *Matrix Computations*,

xtask/src/layers.rs

@@ -37,6 +37,8 @@ impl Class {
 pub const TABLE: &[(&str, Class)] = &[
     // L0: the solver-neutral model.
     ("model", Class::Layer(0)),
+    // L0, no workspace dependencies: deterministic thread pool helpers for the solvers.
+    ("parallel", Class::Standalone),
     // L1: geometry and math kernels.
     ("geometry", Class::Layer(1)),
     ("mesh", Class::Layer(1)),