Module parallel 

The sparse LDU matrix-vector product and its companion vector operations, on the hybrid execution backend — the hot path of every implicit finite-volume solve.

A Krylov solver (conjugate gradient, BiCGStab, GMRES) spends most of its wall clock inside a handful of operations. This module provides each of them as one public entry point that takes a ComputeBackend parameter, in the shape the hybrid-backend epic mandates: dispatch, not a foo_parallel() sibling beside foo().

Operation	Cost	Entry point
sparse product y = A x	O(n_cells + n_faces)	HybridLdu::spmv / HybridLdu::spmv_into
residual r = b - A x	O(n_cells + n_faces)	HybridLdu::residual / HybridLdu::residual_into
scaled residual norm	O(n_cells + n_faces)	HybridLdu::normalised_residual
diagonal reciprocal 1 / diag	O(n_cells)	HybridLdu::diagonal_reciprocal
inner product a . b	O(n)	dot
y := alpha x + y	O(n)	axpy
x := alpha x	O(n)	scale
sqrt(sum x_i^2) / sum abs(x_i)	O(n)	norm_l2 / norm_l1

The solvers that consume them live in crate::krylov (BiCGStab, GMRES) and crate::ldu_matrix::solvers; see crate::krylov::bicgstab_prepared for the entry point that drives this module’s kernels end to end.

The correctness problem this module solves

LduMatrix stores its off-diagonal coefficients per internal face, so the textbook product is a face-based scatter:

for each face f:  y[owner[f]]     += upper[f] * x[neighbour[f]]
                  y[neighbour[f]] += lower[f] * x[owner[f]]

Two faces generally share a cell, so parallelising that face loop directly is a data race: two threads read-modify-write the same y[c]. Where such a thing can be written at all (atomics, unchecked split borrows) it produces silently wrong, run-varying answers.

This module uses the cell-gather reformulation. A one-off index build (LduTopology) inverts the face addressing into a per-cell list of incident faces, turning the product into

for each cell c:  y[c] = diag[c] * x[c]
                       + sum over faces f incident on c of
                             (c is owner ? upper[f] * x[neighbour[f]]
                                         : lower[f] * x[owner[f]])

The loop is now over cells, and every output element y[c] is written by exactly one thread. There is no race, no atomic, and no per-thread scratch buffer.

The two alternatives were considered and rejected. Per-thread partial accumulation costs threads * n_cells extra memory and needs a reduction pass whose association order varies with the schedule, so it is not reproducible. Face colouring needs a graph colouring at build time, still writes each cell once per colour, and its result depends on the colouring produced. The cell-gather index costs one O(n_cells + n_faces) build that is amortised over the thousands of products a solve performs, and it buys exact reproducibility (below).

Determinism

Every kernel in this module returns bit-for-bit identical output on ComputeBackend::Serial and ComputeBackend::CpuMulti, at any thread count, on every run. That is stronger than the usual parallel-reduction contract and it is deliberate, because ComputeBackend::Serial is this workspace’s documented deterministic oracle.

Two separate mechanisms deliver it:

• Products and element-wise kernels (HybridLdu::spmv, HybridLdu::residual, HybridLdu::diagonal_reciprocal, axpy, scale) are bitwise identical also to the pre-existing serial reference LduMatrix::multiply / LduMatrix::residual. LduTopology lists each cell’s incident faces in ascending face index, which is exactly the order in which the serial scatter reaches that cell, so each y[c] accumulates the same additions in the same sequence.
• Reductions (dot, norm_l1, norm_l2, HybridLdu::normalised_residual) sum in fixed-size blocks of REDUCTION_BLOCK elements and then combine the block partials in ascending block order. The association is a function of the array length and REDUCTION_BLOCK alone — never of the thread count or the work-stealing schedule — so a 1-thread run agrees bit for bit with a 64-thread run.

The one thing a reduction here is not is bitwise equal to a flat left-to-right sum such as crate::krylov::vecops::dot or LduMatrix::normalised_residual. Blocked summation reassociates, and floating-point addition is not associative. That difference is small, bounded, and measured rather than asserted — see the “Measured deviation” sections on dot and HybridLdu::normalised_residual. Blocked summation is in fact the more accurate of the two (it is a two-level pairwise-style sum), so this is not a loss of accuracy relative to the flat reference.

When multi-CPU is actually faster

Threading is not free. Below SPMV_MIN_CELLS cells (products) or VECOP_MIN_ELEMENTS elements (vector operations), a CpuMulti request runs the serial kernel on the calling thread instead. Because the two paths are bitwise identical, that size dispatch changes no number a caller can observe — only the wall clock. Both constants were measured on this workspace’s development machine; the tables are on the constants themselves.

For scale, the sparse product on a 262 144-cell / 774 144-face matrix, measured 2026-08-12 on 4 logical cores (spmv_thread_scaling_benchmark, best of 7 samples, time per call):

Worker threads	Time per product	Speed-up
1	6487.26 us	1.00x
2	3481.42 us	1.86x
4	2322.10 us	2.79x
8	1830.46 us	3.54x

The output was asserted bitwise identical at every one of those thread counts, which is the determinism claim above measured rather than argued. Eight workers on four cores still help, which is the signature of a memory-latency-bound gather: extra threads buy more outstanding loads rather than more arithmetic. This is one machine and one mesh; it is not a scaling study, and nothing here has been measured on Android hardware or on a many-core server.

Cargo features

The rayon code paths sit behind the crate’s parallel feature, which is off by default. With the feature off this module still compiles and every entry point still works: ComputeBackend::CpuMulti resolves down to ComputeBackend::Serial via ComputeBackend::resolve and the answer is unchanged. There is no Gpu kernel here yet, so a Gpu request also degrades to the best available CPU path.

Portability

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.

Units

All slices are dimensionless f64 in cell order: element c is the value at cell c, in whatever units the assembled equation carries. No uom typing is applied at this layer, for the same reason crate::krylov::vecops applies none: a Krylov subspace mixes residuals, search directions and solution increments that share no single physical dimension. Units belong on the field and equation layer that assembles the matrix, and are not stripped there.

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;

// 3-cell symmetric tridiagonal  [[2,-1,0],[-1,2,-1],[0,-1,2]]
let mut m = LduMatrix::new(3, vec![0, 1], vec![1, 2]);
m.diag = vec![2.0, 2.0, 2.0];
m.upper = vec![-1.0, -1.0];
m.lower = vec![-1.0, -1.0];
let m = Arc::new(m);

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

assert_eq!(ldu.spmv(&x, ComputeBackend::Serial), vec![1.0, 0.0, 1.0]);

// Asking for multi-CPU gives a bit-for-bit identical answer, whether or not
// the `parallel` feature is compiled in.
assert_eq!(
    ldu.spmv(&x, ComputeBackend::CpuMulti),
    ldu.spmv(&x, ComputeBackend::Serial),
);

Structs

HybridLdu

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

LduTopology

The face addressing of an LduMatrix, inverted into per-cell incident-face lists so the matrix-vector product can be parallelised over cells.

Constants

CELL_BLOCK

Number of cells in one block of the cell-parallel kernels.

REDUCTION_BLOCK

Number of elements summed serially inside one block of a reduction, before block partials are combined.

SPMV_MIN_CELLS

Cell count below which a ComputeBackend::CpuMulti request runs the serial product kernel instead.

VECOP_MIN_ELEMENTS

Element count below which a ComputeBackend::CpuMulti request runs the serial kernel instead, for the vector operations dot, axpy, norm_l1 and norm_l2.

Functions

axpy

AXPY update y := alpha * x + y, in place, on the chosen backend.

dot

Inner product sum_i a_i * b_i, on the chosen backend.

norm_l1

L1 norm sum_i abs(x_i), on the chosen backend.

norm_l2

L2 (Euclidean) norm sqrt(sum_i x_i^2), on the chosen backend.

scale

Scale x := alpha * x, in place, on the chosen backend.

spmv_backend_for

The ComputeBackend this module would actually use for a sparse product over n_cells cells if asked for requested — without running anything.

vecop_backend_for

The ComputeBackend this module would actually use for a vector operation over n elements if asked for requested — without running anything.