Struct HybridLdu 

pub struct HybridLdu { /* private fields */ }

An LduMatrix bundled with its cell-gather index, exposing the per-iteration kernels a Krylov solver needs on any ComputeBackend.

Construct once per assembled matrix and call the kernels many times: the O(n_cells + n_faces) index build happens in Self::new, and every kernel call after that is index-free. Both the matrix and the index are held behind Arc, so cloning is cheap and the value can be shared across threads.

Which backend runs

Every kernel takes the backend as a parameter — there is no _parallel sibling API. What actually runs is spmv_backend_for / vecop_backend_for applied to the request: an unavailable backend degrades, Gpu degrades (no GPU kernel here yet), and a problem below the measured size floor runs serially. None of those degradations changes the answer. Self::auto_backend applies the crate-wide policy crate::compute::select_backend if you would rather not choose.

Thread pool

The parallel kernels use rayon’s ambient pool — the global one by default. No pool is built here, and none is built per call. A caller that wants a specific worker count builds its own rayon::ThreadPool and calls these kernels inside its install(...) scope; the parallel iterators then run on that pool. Because every kernel is bitwise deterministic at any thread count, that choice affects wall time only.

Units

All vectors are dimensionless f64 in cell order; see the module documentation.

Example

use std::sync::Arc;
use outram_foam_basic_lib::compute::ComputeBackend;
use outram_foam_basic_lib::ldu_matrix::LduMatrix;
use outram_foam_basic_lib::ldu_matrix::parallel::HybridLdu;

let mut m = LduMatrix::new(4, vec![0, 1, 2], vec![1, 2, 3]);
m.diag = vec![4.0, 4.0, 4.0, 4.0];
m.upper = vec![-1.0, -1.0, -1.0];
m.lower = vec![-1.0, -1.0, -1.0];
let m = Arc::new(m);

let ldu = HybridLdu::new(Arc::clone(&m));
let x = vec![1.0, 2.0, 3.0, 4.0];

// Bit-for-bit agreement with the pre-existing serial reference kernel.
assert_eq!(ldu.spmv(&x, ComputeBackend::CpuMulti), m.multiply(&x));

// Reassembling the same mesh reuses the index instead of rebuilding it.
let mut m2 = (*m).clone();
m2.diag = vec![8.0, 8.0, 8.0, 8.0];
let ldu2 = ldu.with_matrix(Arc::new(m2)).expect("same addressing");
assert_eq!(ldu2.matrix().diag[0], 8.0);

Implementations

impl HybridLdu

pub fn new(matrix: Arc<LduMatrix>) -> Self

Build the cell-gather index for matrix and wrap it for hybrid execution.

The O(n_cells + n_internal_faces) index build happens here, once.

Arguments

• matrix — the assembled sparse system, shared by Arc so the caller may keep its own handle.

Panics

Panics if the matrix addressing is inconsistent — see LduTopology::from_matrix.

pub fn with_matrix(&self, matrix: Arc<LduMatrix>) -> Option<Self>

Reuse this index with a reassembled matrix over the same mesh.

A finite-volume solver rebuilds coefficients every outer iteration while the mesh addressing is unchanged. This returns a HybridLdu over the new coefficients that shares the existing index — an O(n_cells + n_internal_faces) addressing check instead of a full rebuild.

Returns

None if matrix does not have exactly the same owner/neighbour addressing (see LduTopology::matches), in which case the caller must go through Self::new. Returning None rather than reusing a mismatched index is deliberate: a stale index would produce a wrong answer, not a crash.

pub fn matrix(&self) -> &Arc<LduMatrix>

The matrix these kernels operate on.

pub fn topology(&self) -> &Arc<LduTopology>

The cell-gather index, shared behind an Arc.

pub fn auto_backend(&self) -> ComputeBackend

The backend the crate-wide policy crate::compute::select_backend picks for a product on this matrix.

Passes n_cells as the work-item count, because the cell-gather kernel’s independent work items are cells. The returned backend is guaranteed available, so it can be handed straight to Self::spmv.

Note this is the crate-wide policy, whose CpuMulti threshold is documented as a placeholder; this module’s own measured threshold SPMV_MIN_CELLS is applied afterwards by the kernel itself, so a too-eager CpuMulti from here still runs serially below that size.

pub fn spmv_into(&self, x: &[f64], y: &mut [f64], backend: ComputeBackend)

Sparse matrix-vector product y = A x, writing into a caller-owned buffer.

This is the kernel a Krylov solver calls once (BiCGStab: twice) per iteration and where it spends most of its time. Prefer it over Self::spmv inside a solver loop — it reuses y instead of allocating a fresh Vec every iteration.

Arguments

• x — input vector in cell order, exactly n_cells long, dimensionless.
• y — output buffer, exactly n_cells long. Fully overwritten; its prior contents are irrelevant. Aliasing with x is impossible because y is a unique borrow.
• backend — requested execution backend; see spmv_backend_for for what will actually run.

Determinism

Bitwise identical to LduMatrix::multiply on every backend and at any thread count. See the module documentation, “Determinism”.

Panics

Panics if x.len() or y.len() differs from the matrix’s n_cells.

pub fn spmv(&self, x: &[f64], backend: ComputeBackend) -> Vec<f64>

Sparse matrix-vector product y = A x, allocating the result.

Convenience wrapper over Self::spmv_into with identical semantics and the same bitwise guarantee. Inside a solver loop prefer Self::spmv_into, which does not allocate.

Arguments

• x — input vector in cell order, exactly n_cells long.
• backend — requested execution backend.

Returns

A fresh Vec<f64> of length n_cells.

Panics

Panics if x.len() differs from the matrix’s n_cells.

pub fn residual_into( &self, x: &[f64], b: &[f64], r: &mut [f64], backend: ComputeBackend, )

Residual r = b - A x, writing into a caller-owned buffer.

Arguments

• x — current solution estimate, n_cells long.
• b — the equation’s source vector, n_cells long.
• r — output buffer, n_cells long. Fully overwritten.
• backend — requested execution backend.

Determinism

Bitwise identical to LduMatrix::residual on every backend and at any thread count: the product is bitwise exact and the subtraction is element-wise, so no reassociation can occur.

Panics

Panics if any of x, b, r has a length other than n_cells.

pub fn residual( &self, x: &[f64], b: &[f64], backend: ComputeBackend, ) -> Vec<f64>

Residual r = b - A x, allocating the result.

Convenience wrapper over Self::residual_into; identical semantics and the same bitwise guarantee against LduMatrix::residual.

Panics

Panics if x.len() or b.len() differs from n_cells.

pub fn normalised_residual( &self, x: &[f64], b: &[f64], backend: ComputeBackend, ) -> f64

OpenFOAM-style scaled residual norm, ||b - A x||_1 / (0.5 * (||A x||_1 + ||b||_1)).

This is the convergence measure the solvers in crate::ldu_matrix::solvers test against a tolerance, so it is evaluated once per iteration. The scaling makes the number dimensionless and roughly comparable across equations of very different magnitude. When the denominator underflows below f64::EPSILON the unscaled ||b - A x||_1 is returned instead, matching LduMatrix::normalised_residual.

Arguments

• x — current solution estimate, n_cells long.
• b — source vector, n_cells long.
• backend — requested execution backend.

Returns

A dimensionless non-negative scalar.

Determinism, and how this differs from the serial reference

This is a reduction, so unlike the product it cannot be bitwise identical to a flat left-to-right sum. It sums in blocks of REDUCTION_BLOCK elements and combines the block partials in ascending block order, an association fixed by n_cells alone. Therefore repeated runs agree bit for bit, and a 1-thread run agrees bit for bit with an N-thread run — but the value differs from LduMatrix::normalised_residual in the last bits.

Measured deviation. Methodology: a fixed-seed xorshift64* pseudorandom, diagonally dominant symmetric 7-point-stencil matrix on a 32x32x32 mesh (32 768 cells, 95 232 internal faces), with pseudorandom x and b uniform on [-1, 1); compare this function against LduMatrix::normalised_residual on the same inputs. Pass criterion: relative difference <= 1e-13, plus bitwise equality between the two backends.

Result, measured 2026-08-12 by the test normalised_residual_matches_flat_reference in parallel/tests.rs: blocked 1.74107987267472608e0, flat 1.74107987267473940e0, relative difference 7.6520e-15 — about 34 units in the last place, an order of magnitude inside the gate — and bitwise equality across backends held.

Interpretation: the convergence measure a solver tests against a tolerance is unchanged for any practical tolerance. Solver tolerances in this crate are 1e-6 to 1e-12 on a quantity whose two summation orders differ in the fifteenth digit, so switching backend cannot change an iteration count.

Panics

Panics if x.len() or b.len() differs from n_cells.

pub fn diagonal_reciprocal(&self, backend: ComputeBackend) -> Vec<f64>

Element-wise reciprocal of the diagonal, 1 / diag[c] for every cell.

This is the Jacobi preconditioner, and the first step of the DIC/DILU preconditioner setups. It is recomputed whenever the matrix is reassembled, which for a transient solve is every timestep. O(n_cells).

A zero diagonal yields an infinity and a NaN diagonal propagates, exactly as the scalar expression 1.0 / diag[c] would. No clamping or masking is applied, because silently substituting a finite value would hide a singular assembly. Callers needing a guard should test the result with f64::is_finite.

Arguments

• backend — requested execution backend.

Returns

A fresh Vec<f64> of length n_cells, in cell order. Its units are the reciprocal of the diagonal’s units.

Determinism

Bitwise identical on every backend and at any thread count: each element is one independent division.

Example

use std::sync::Arc;
use outram_foam_basic_lib::compute::ComputeBackend;
use outram_foam_basic_lib::ldu_matrix::LduMatrix;
use outram_foam_basic_lib::ldu_matrix::parallel::HybridLdu;

let mut m = LduMatrix::new(3, vec![], vec![]);
m.diag = vec![2.0, 4.0, 0.0];
let ldu = HybridLdu::new(Arc::new(m));

let rd = ldu.diagonal_reciprocal(ComputeBackend::CpuMulti);
assert_eq!(rd[0], 0.5);
assert_eq!(rd[1], 0.25);
assert!(rd[2].is_infinite()); // a singular row is reported, not hidden

pub fn diagonal_reciprocal_into(&self, out: &mut [f64], backend: ComputeBackend)

Element-wise reciprocal of the diagonal, into a caller-owned buffer.

Same semantics as Self::diagonal_reciprocal without the allocation.

Arguments

• out — output buffer, exactly n_cells long. Fully overwritten.
• backend — requested execution backend.

Panics

Panics if out.len() differs from n_cells.

Trait Implementations

impl Clone for HybridLdu

fn clone(&self) -> Self

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Debug for HybridLdu

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.