Module parallel 

Batched scalar root finding on the hybrid execution backend — solve N independent f(x) = 0 problems at once, serially or across CPU cores.

The shape of the problem

A finite-volume solver almost never wants one root. It wants one root per cell: given a cell-wise target enthalpy h_target[c], find the temperature T[c] satisfying h(T[c]) - h_target[c] = 0, for every cell in the mesh. Those N problems are completely independent of one another — no cell’s iteration reads another cell’s iterate — which is what makes the batch worth threading at all.

This module therefore provides batched entry points only. A single scalar root solve should stay serial and inline in its caller: one root is a few microseconds of work and threading it is a straight loss.

Operation	Derivative needed?	Entry point
bisection or Brent over a bracket	no	solve_bracketed_batch
safeguarded Newton over a bracket	yes	solve_newton_batch
closed-form roots of a x + b	—	linear_roots_batch
closed-form roots of a x^2 + b x + c	—	quadratic_roots_batch
closed-form roots of a x^3 + b x^2 + c x + d	—	cubic_roots_batch

Hybrid means dispatch, not two APIs

Every entry point takes a ComputeBackend parameter. There is no solve_batch_parallel() sibling beside solve_batch(). With the parallel feature off, ComputeBackend::CpuMulti resolves down to ComputeBackend::Serial through ComputeBackend::resolve and the answer is unchanged — bit for bit, not merely close. There is no Gpu kernel here yet, so a Gpu request degrades to the best available CPU path.

Determinism — bitwise identical, and why that is stronger here

Every kernel in this module returns bit-for-bit identical output on ComputeBackend::Serial and ComputeBackend::CpuMulti, at any thread count, on every run — provided the caller’s residual closure is itself a deterministic pure function of (index, x).

This is easier to achieve than it is for a reduction and the reason is worth stating, because it is the thing most people get wrong about parallel root finding. A parallel sum has to re-associate a + b + c, and floating-point addition is not associative, so its last bits depend on the split. A batch of root solves has no cross-element arithmetic at all: lane i’s answer is a pure function of lane i’s bracket and closure evaluations. Splitting the lanes across threads cannot perturb any of them. Both backends call the very same #[inline] per-lane kernel; only the identity of the calling thread differs.

The one way a caller can break this is to supply a closure that is not a pure function — one that reads a random number generator, or accumulates into shared interior-mutable state, or depends on the calling thread. Such a closure makes the batch non-reproducible no matter what this module does. The Sync bound permits it; the documented contract forbids it.

Verified by the bitwise_* tests in parallel/tests.rs, which compare serial against rayon pools of 1, 2, 4 and 8 workers on a batch deliberately built with wildly varying per-lane iteration counts.

Load imbalance — why there is no hand-rolled partition

Iteration counts vary per lane, sometimes by an order of magnitude: a lane whose bracket happens to straddle the root tightly converges in three iterations, and its neighbour takes eighty. A static equal-split across P threads therefore ends up waiting on whichever chunk drew the hard lanes.

The parallel paths here use rayon’s adaptive splitting (par_iter_mut, into_par_iter) with no min_len floor on the iterative solvers, so an idle worker steals down to individual lanes when the work is skewed. That is a deliberate choice, not an oversight: the whole reason to prefer work-stealing here is the per-lane cost variance. The closed-form polynomial kernels are the exception — every lane there costs the same handful of floating-point operations, so they set a granularity floor of POLY_BLOCK lanes to stop the splitter subdividing below a task that cannot pay for itself. That floor affects wall clock only; it cannot change a value, because each lane is independent.

For scale, a 65 536-lane deliberately-imbalanced batch (half the lanes given a 2e-6-wide bracket, half given [0, 1e6]) under Brent, measured 2026-08-12 on 4 logical cores by the #[ignore]d root_batch_thread_scaling_benchmark, best of 7 samples, with a second independent run alongside:

Worker threads	Time	Speed-up	(repeat)	Bitwise vs serial
serial reference	19322.52 us	1.00x	1.00x	—
1	20510.52 us	0.94x	0.96x	identical
2	10123.65 us	1.91x	1.88x	identical
4	7075.89 us	2.73x	3.66x	identical
8	7314.01 us	2.64x	3.57x	identical

The “identical” column is the determinism claim above measured rather than argued, and it is asserted by the benchmark itself, not merely printed. The one-worker row costs about 6% against plain serial, which is the price of going through rayon at all. Unlike the memory-latency-bound sparse product in crate::ldu_matrix::parallel, eight workers on four cores buy nothing further here — the expected signature of a compute-bound kernel that already saturates the cores it has. This is one machine, one batch and two runs; it is not a scaling study, and nothing here has been measured on Android hardware or on a many-core server.

Non-convergence is reported, never swallowed

A batch of 10 000 problems in which 3 fail must say so. This module reports failure per lane, and makes the failure impossible to ignore by construction:

• Every lane carries a RootStatus and is reachable through RootBatch::solutions.
• RootSolution::root returns Option<f64> — Some(x) only when that lane converged. There is no accessor that hands out a bare f64 without the caller having seen the status; the diagnostic iterate is behind the deliberately-named RootSolution::last_iterate.
• RootBatch::roots is the all-or-nothing path: it returns Err(RootBatchFailure) naming the failure count and the first failing lane, rather than a plausible-looking Vec<f64>.

A non-converged lane is never clamped to a bracket endpoint and called a root. A lane that ran out of iterations reports its best iterate with RootStatus::MaxIterations; a lane whose bracket does not straddle a root reports RootStatus::NotBracketed and a NaN iterate, because there is no honest number to return.

Why this module lives under math/ and not polynomial/

The bulk of it is a general scalar root finder over a caller-supplied residual — h(T) - h_target, a Peng-Robinson compressibility, a saturation curve. None of that is polynomial. math/ is already the home of the crate’s iterative inversions ([crate::math::inv_inc_gamma] is a Newton-refined inverse), so a general solver belongs beside them; a caller inverting an enthalpy would not think to look in a module called polynomial. The batched closed-form polynomial kernels are thin wrappers over crate::polynomial’s existing per-equation solvers and are kept here with the rest of the batching so that the crate has one batched root-finding vocabulary rather than two dialects in two modules.

Units

Everything here is dimensionless f64, for the same reason crate::ldu_matrix::parallel and crate::krylov::vecops are: a generic root finder has no single physical dimension. The abscissa x of one lane may be a temperature in kelvin and of another a compressibility factor, and the residual’s dimension is whatever the caller’s function returns.

uom typing is not stripped to get here — it is applied at the boundary, by the caller, exactly as the hybrid-backend epic requires: convert into the batch, convert back out. The doctest on solve_newton_batch shows that boundary explicitly, inverting a uom-typed enthalpy relation for a ThermodynamicTemperature.

Cargo features and portability

The rayon paths sit behind the crate’s parallel feature, which is off by default; with it off this module still compiles and every entry point still works. rayon is pure Rust with no system component, so everything here compiles and runs on aarch64-linux-android / Termux exactly as on desktop. Nothing in this module is target-gated.

Example

use outram_foam_basic_lib::compute::ComputeBackend;
use outram_foam_basic_lib::math::parallel::{
    solve_bracketed_batch, RootMethod, RootProblem, RootSettings,
};

// 4 independent problems: x^2 = k for k = 1, 2, 3, 4, each bracketed on [0, 4].
let targets = [1.0, 2.0, 3.0, 4.0];
let problems: Vec<RootProblem> = (0..4).map(|_| RootProblem::new(0.0, 4.0)).collect();

let batch = solve_bracketed_batch(
    &problems,
    RootMethod::Brent,
    RootSettings::default(),
    ComputeBackend::CpuMulti,
    |i, x| x * x - targets[i],
);

let roots = batch.roots().expect("all four lanes converge");
assert!((roots[1] - 2.0_f64.sqrt()).abs() < 1e-12);

// Asking for multi-CPU gives a bit-for-bit identical answer, whether or not
// the `parallel` feature is compiled in.
let serial = solve_bracketed_batch(
    &problems,
    RootMethod::Brent,
    RootSettings::default(),
    ComputeBackend::Serial,
    |i, x| x * x - targets[i],
);
assert_eq!(roots, serial.roots().unwrap());

Structs

RootBatch

A batch of N lane outcomes, in the same order as the problems handed in.

RootBatchFailure

One or more lanes of a RootBatch did not converge.

RootProblem

One bracketed root problem: an interval believed to contain a root, plus an optional starting iterate for the Newton solver.

RootSettings

Convergence tolerances and the iteration cap, shared by both iterative solvers.

RootSolution

The outcome of a single lane: its status, its iterate, its residual and the work it took.

Enums

RootMethod

Which derivative-free bracketed method solve_bracketed_batch should run.

RootStatus

How one lane of a batch ended.

Constants

POLY_BLOCK

Number of lanes in one task of the closed-form polynomial kernels (linear_roots_batch, quadratic_roots_batch, cubic_roots_batch).

POLY_ROOTS_MIN_EQUATIONS

Equation count below which a ComputeBackend::CpuMulti request runs the closed-form polynomial kernels on the calling thread instead.

ROOT_BATCH_MIN_PROBLEMS

Problem count below which a ComputeBackend::CpuMulti request runs the iterative root-finding kernels (solve_bracketed_batch, solve_newton_batch) on the calling thread instead.

Functions

cubic_roots_batch

Roots of N independent cubic equations a x^3 + b x^2 + c x + d = 0, on the chosen backend.

linear_roots_batch

Roots of N independent linear equations a x + b = 0, on the chosen backend.

poly_roots_backend_for

The ComputeBackend the closed-form polynomial kernels would actually use for n equations if asked for requested — without running anything.

quadratic_roots_batch

Roots of N independent quadratic equations a x^2 + b x + c = 0, on the chosen backend.

root_batch_backend_for

The ComputeBackend the iterative solvers would actually use for n problems if asked for requested — without running anything.

solve_bracketed_batch

Solve N independent bracketed root problems without a derivative, on the chosen backend.

solve_newton_batch

Solve N independent root problems by bracket-safeguarded Newton iteration, on the chosen backend.