HyBIT documentation

HyBIT 0.8.0 is the current published release. The immutable release source is tagged v0.8.0. Rendered documentation: https://michioga.github.io/hybit/

The site is built with mdBook from the Markdown files in this directory and GitHub Pages deploys the frozen release documentation from main.

Start here

For library users:

1. User guide — choose a problem class, configure a solver, prepare reusable state, solve, and validate the result.
2. Solver selection and Krylov theory — understand why SPD uses PCG, GeneralSquare uses FGMRES, how preconditioning changes the mathematical contract, and what MINRES/LSQR-class future routes require.
3. Prepared execution — understand analyze -> prepare -> solve-many, reuse boundaries, and memory ownership.
4. General-square systems — FGMRES, restart policies, Jacobi versus ILU(0), pivot stabilization, and current limitations.
5. Structural Auto API — 3-D structural SPD systems with rigid-body coarse correction.
6. Threading and host-runtime interoperability — Rayon/OpenMP thread-count precedence, MPI/OpenMP embedding, and oversubscription rules.

For numerical background:

• Solver selection and Krylov theory — matrix-class contracts, PCG/FGMRES mechanics, restart cost, compatibility, and breakdown.
• Hybrid SPD mathematics — selective local direct correction, overlap weighting, and resumable PCG logic.
• Architecture — crate boundaries, execution architecture, resident CPU work, FGMRES checkpoints, and internal policy flow.
• ABTM topology algebra — planned symbolic/topology layer, metadata-first numerical pruning, adaptive/block layouts, and validation sequence for CPU, GPU, and distributed execution.

For project status:

• Development status — published release versus the current development branch.
• Roadmap — completed and planned work.
• FEM benchmark and benchmark README — reproducible validation.
• Publishing — release discipline and immutable release tags.

Documentation layers

Benchmark results are evidence for tested cases, not universal performance claims.

HyBIT user guide

Development documentation for develop/0.8.0.

HyBIT 0.7.0 remains the current published crates.io release. The 0.8 APIs described here are development APIs and may change before the next release.

1. Choose the mathematical problem class

HyBIT does not infer the mathematical class of an arbitrary matrix. The caller declares the contract with MatrixProblemClass.

For the mathematical reason behind these routes, including PCG positivity, FGMRES flexible right preconditioning, restart memory/work tradeoffs, and the future MINRES contract, see Solver selection and Krylov theory. Spd is the default. Do not label a matrix SPD merely to reach PCG: the algorithm and SPD preconditioners rely on symmetry and positive definiteness.

Rectangular systems are not currently routed by the high-level solver.

CSR32 input contract

Csr32Matrix::new validates the storage contract:

• row_ptr.len() == nrows + 1;
• row_ptr[0] == 0;
• row_ptr is monotonically nondecreasing;
• row_ptr[nrows] == nnz;
• col_idx.len() == values.len();
• every column index is smaller than ncols;
• all matrix values are finite;
• nnz <= u32::MAX.

CSR rows do not have to be column-sorted and duplicate column entries are allowed. Matrix-vector multiplication therefore follows the caller’s stored row order. diagonal() sums duplicate diagonal entries.

Problem-specific requirements are stricter than the storage contract. For example, the SPD route requires a complete positive diagonal. GeneralSquare Jacobi requires a complete finite nonzero diagonal. GeneralSquare ILU(0) requires diagonal entries to exist, then privately canonicalizes rows before factorization.

Defaults at a glance

The most important high-level defaults on the 0.8 development branch are:

Defaults are conservative compatibility choices, not claims that they are optimal for every matrix family.

2. Use the smallest correct interface

One-shot SPD:

use hybit::{Csr32Matrix, HybitSolver};

let a = Csr32Matrix::new(
    3,
    3,
    vec![0, 2, 5, 7],
    vec![0, 1, 0, 1, 2, 1, 2],
    vec![2.0, -1.0, -1.0, 2.0, -1.0, -1.0, 2.0],
)?;
let b = vec![1.0, 0.0, 1.0];
let mut x = vec![0.0; 3];

let solver = HybitSolver::new();
let report = solver.solve_csr32(&a, &b, &mut x)?;
assert!(report.converged());

GeneralSquare with opt-in ILU(0):

use hybit::{
    GeneralSquareOptions, GeneralSquarePreconditionerPolicy, HybitSolver,
    MatrixProblemClass, SolverOptions,
};

let mut solver = HybitSolver::new();
solver.set_problem_class(MatrixProblemClass::GeneralSquare);
solver.set_options(SolverOptions {
    relative_tolerance: 1.0e-8,
    absolute_tolerance: 0.0,
    max_iterations: 1800,
})?;

solver.set_general_square_preconditioner_policy(
    GeneralSquarePreconditionerPolicy::Ilu0,
);

// Measured E5 setting, not an automatic default.
solver.set_general_square_options(GeneralSquareOptions {
    restart: 3,
    ..GeneralSquareOptions::default()
})?;

let mut x = vec![0.0; a.nrows()];
let report = solver.solve_csr32(&a, &b, &mut x)?;

Run the complete example with:

cargo run --release -p hybit --example general_square

Configure before analyze

For a reproducible lifecycle, set the problem class, solver tolerances, backend/execution policy, GeneralSquare preconditioner, and restart policy before calling analyze_csr32.

The analysis records the problem class, execution target/policy, backend choice, and matrix signatures. prepare_csr32 validates compatible analysis state and then captures the current solver/preconditioner options into the prepared context. Changing the HybitSolver object afterwards does not mutate an already-created prepared context.

3. Prefer prepared execution for repeated right-hand sides

let analysis = solver.analyze_csr32(&a)?;
let mut prepared = solver.prepare_csr32(&a, &analysis)?;

let mut x1 = vec![0.0; a.nrows()];
let first = prepared.solve(&a, &b1, &mut x1)?;

let mut x2 = vec![0.0; a.nrows()];
let second = prepared.solve(&a, &b2, &mut x2)?;

Lifecycle:

matrix + solver policy
        |
        v
     analyze
        |
        |  validate problem contract
        |  profile matrix
        |  resolve backend/execution policy
        |  record structure/value signatures
        v
     prepare
        |
        |  build matrix-dependent preconditioner state
        |  allocate reusable Krylov workspace
        v
     solve #1
        |
        v
     solve #2 ... N
        |
        +-- reuse matrix-dependent state

Prepared contexts currently require unchanged CSR structure and unchanged coefficient bit patterns. See PREPARED_API.md.

4. Validate the result

At minimum, inspect:

• report.converged();
• report.relative_residual;
• report.iterations;
• report.solver;
• report.preconditioner;
• report.backend.

report.final_residual is the absolute Euclidean residual norm. For nonzero RHS,

report.relative_residual = ||b - A x||_2 / ||b||_2

For a zero RHS, relative_residual is reported as the absolute final residual instead of dividing by zero.

For engineering/scientific work, independently recompute the same residual from the original matrix and RHS. Linear-solver convergence does not validate the physical model, mesh, units, boundary conditions, or material law.

5. GeneralSquare choices are independent

Preconditioner:

• Jacobi — default, low setup/storage, finite nonzero diagonal required.
• Ilu0 — opt-in, stronger coupling model, canonical CSR factorization, currently serial triangular application.

Restart:

• Fixed — default; restart 30.
• Escalating — bounded staged restart growth.
• BudgetAware — restart-boundary residual controller.

E5 found FGMRES(3)+ILU(0) fastest among tested restart values on three synthetic hard families. HyBIT does not make that combination the global default because ordering and matrix family can materially change ILU(0) quality.

See GENERAL_SQUARE.md.

6. SPD and structural paths

Generic SPD uses solve_csr32 / prepare_csr32 and PCG-based policies.

For 3-D structural SPD systems with node coordinates, use solve_structural_csr32 / prepare_structural_csr32 to enable the rigid-body coarse path. See STRUCTURAL_AUTO_API.md.

7. Backend and execution scope

CSR32 is the public sparse input representation. ABTM is an internal execution/topology representation.

GeneralSquare FGMRES currently accepts ExecutionPolicy::Auto or ExecutionPolicy::Cpu. Resident Rayon and GPU routing are not active for this path yet.

8. Rust versus C/C++/Fortran

The repository retains the stable 0.7 C ABI and language examples. The 0.8 Rust API can select GeneralSquare, FGMRES restart policy, and ILU(0).

The current C ABI can report FGMRES/ILU(0) enum codes, but does not yet expose setters for selecting the GeneralSquare problem class and its preconditioner. Treat the new GeneralSquare path as Rust-API-only until that surface is added and validated.

9. Evaluation sequence for a new matrix family

1. confirm the mathematical problem class;
2. run a simple baseline;
3. independently verify the true residual;
4. record iterations, setup, solve time, and persistent state;
5. change one policy at a time;
6. for ILU(0), test ordering sensitivity;
7. test solve-many reuse separately from warm starts;
8. keep conclusions scoped to tested matrices and hardware.

Solver selection and Krylov theory

Development documentation for develop/0.8.0.

This document describes the mathematical contracts behind HyBIT’s current square-system routing. It distinguishes implemented paths from future routes.

1. Start from the matrix class, not the solver name

For a linear system

A x = b

the first solver decision is a statement about A, not a performance tuning choice.

                         square system
                              |
             +----------------+----------------+
             |                                 |
        symmetric?                             no / unknown
             |                                 |
        +----+----+                            v
        |         |                    GeneralSquare
   positive     indefinite                    |
   definite        |                          v
        |          |                       FGMRES
        v          v                    + Jacobi / ILU(0)
       SPD      SymmetricIndefinite
        |          |
        v          v
       PCG       MINRES [future]
        |
        +-- Jacobi
        +-- Hybrid local direct
        +-- algebraic two-level
        +-- structural rigid-body two-level

Current high-level routing is:

HyBIT does not infer positive definiteness from a positive diagonal, symmetry flag, or successful factorization. MatrixProblemClass is a caller-declared mathematical contract.

A singular, unconstrained, or nearly singular system is not made suitable for a method merely by choosing one of these labels. For example, an unconstrained structural stiffness matrix with rigid-body null modes is not SPD.

2. Why SPD systems use PCG

For PCG, the intended operator contract is

A = A^T
x^T A x > 0   for every nonzero x

and the preconditioner must preserve the corresponding positive-definite inner-product structure.

With an SPD preconditioner M, preconditioned CG can be interpreted through an SPD transformed problem such as

M^(-1/2) A M^(-1/2) y = M^(-1/2) b
x = M^(-1/2) y

in exact arithmetic. This structure is what supports the short recurrence and the conjugacy/orthogonality relations on which CG depends.

The important consequence is operational:

A preconditioner is not interchangeable merely because it approximately solves A z = r. PCG also needs the preconditioned recurrence to retain the required symmetry and positivity.

HyBIT’s PCG implementation checks two quantities that directly expose loss of those assumptions:

r^T M^-1 r > 0
p^T A p    > 0

A non-finite or non-positive value produces numerical breakdown rather than silently continuing a CG recurrence whose mathematical basis has failed.

PCG workspace and per-iteration shape

The established PcgWorkspace owns five full-size f64 vectors:

A*x / scratch
r
z = M^-1 r
p
A*p

so its numeric payload is

5 * n * sizeof(f64) = 40 n bytes

before counting matrix and preconditioner storage.

A normal PCG iteration has one operator application, one preconditioner application, a small number of global reductions, and vector updates. Its Krylov storage therefore stays O(n) rather than growing with iteration count.

3. Why HyBIT may continue or restart PCG

HyBIT’s SPD controller divides a solve into observation and strengthening stages. There are two mathematically different cases.

Same preconditioner: continue the recurrence

If a probe or telemetry boundary does not change A or M, the PcgSession retains the current search direction and r^T M^-1 r state.

PCG iterations
      |
      +---- telemetry boundary ----+
      |                            |
      |       A and M unchanged    |
      +---------------------------> continue same recurrence

This is deliberately different from calling PCG again with the current x, which would discard Krylov information.

Changed preconditioner: restart PCG

If selective-direct Hybrid strengthening changes the preconditioner, the old CG conjugacy relation no longer applies to the new operator/preconditioner pair. HyBIT therefore starts a new PCG recurrence from the current solution.

PCG with M0
    |
    +-- difficult region detected
    |
    v
build stronger M1
    |
    v
restart PCG from current x

This restart rule is a mathematical requirement of the current PCG design, not just a controller implementation detail.

See HYBRID_MATH.md for the SPD construction used by the Hybrid preconditioner.

4. SPD preconditioners in the current routes

The generic SPD and structural paths use preconditioners designed to preserve PCG-compatible structure.

Jacobi

The SPD Jacobi path requires a complete positive diagonal. Its approximation is weak but cheap and positive definite when those requirements are satisfied.

Hybrid local direct

For an SPD principal submatrix

A_k = A[H_k, H_k]

HyBIT uses local Cholesky factors and symmetric overlap weighting. The local terms are assembled in a form intended to remain symmetric positive semidefinite, while uncovered DOFs retain positive Jacobi contribution.

Algebraic two-level

The optional generic algebraic coarse correction is an SPD/PCG path. It is disabled by default because an arbitrary generic matrix does not necessarily carry the blocked/geometry semantics of a structural model.

Structural rigid-body two-level

The explicit structural route combines 3x3 block-Jacobi with a six-rigid-body- mode coarse space. It remains an SPD/PCG route; geometry is used to construct the coarse correction, not to change the Krylov method.

5. Why GeneralSquare uses FGMRES

A general real square matrix may be nonsymmetric, nonnormal, or indefinite. The SPD inner-product assumptions behind PCG therefore cannot be presumed.

HyBIT routes MatrixProblemClass::GeneralSquare through restarted, right-preconditioned FGMRES.

At Arnoldi step j:

v_j                  orthonormal Arnoldi basis
z_j = M_j^-1 v_j     independently stored preconditioned basis
w   = A z_j

After orthogonalization, the flexible Arnoldi relation is conceptually

A Z_m = V_(m+1) Hbar_m

and the update is formed from Z, not directly from V:

x_m = x_0 + Z_m y_m

This separation is the key distinction for flexible preconditioning. If M_j changes with iteration, z_j cannot in general be reconstructed later from v_j using one fixed inverse.

The current high-level GeneralSquare route still prepares a fixed Jacobi or fixed ILU(0). FGMRES nevertheless keeps V and Z distinct, so the low-level recurrence already supports a future iteration-dependent preconditioner through FlexiblePreconditioner.

See GENERAL_SQUARE.md for the current Jacobi/ILU(0) preconditioner policies.

6. Restart is a memory/work/convergence tradeoff

Unrestarted GMRES-style methods retain an expanding Krylov basis. Restart limits that state to a maximum Arnoldi dimension m, but discards the old basis at a cycle boundary after updating x.

For problem size n and restart m, the current FgmresWorkspace owns:

V        (m + 1) * n
Z              m * n
r, w           2 * n
H        (m + 1) * m
cs, sn          2 * m
g              (m + 1)
y                    m

Therefore FgmresWorkspace::bytes() counts the numeric payload

8 * [ (2m + 3)n + m^2 + 5m + 1 ] bytes

for f64, excluding Rust Vec headers and allocator overhead.

The basis part is O(nm). Because the implementation uses two-pass modified Gram-Schmidt, orthogonalization work over a full restart cycle grows roughly as O(n m^2).

Consequently:

• larger m retains a richer Krylov subspace and can reduce restart damage;
• larger m consumes more memory;
• larger m increases orthogonalization work;
• a stronger preconditioner can make a small restart faster even when it takes more iterations.

This is why E5’s observed FGMRES(3)+ILU(0) result is treated as matrix-family evidence, not as a universal default.

7. What the restart policies change

The GeneralSquare preconditioner choice and restart policy are independent.

Fixed

Uses one restart dimension for every cycle. The current default is 30.

Escalating

Starts from the configured restart and increases it in bounded stages up to max_restart.

BudgetAware

Observes exact residuals at restart boundaries and can increase the next cycle’s restart dimension when the remaining iteration budget is under sustained pressure.

All policies reuse one prepared FgmresWorkspace whose capacity is large enough for the maximum requested restart. Changing the next cycle’s restart does not require reallocating the Arnoldi basis.

8. True residual versus Hessenberg estimate

Inside a GMRES cycle, Givens rotations provide a cheap residual estimate from the small Hessenberg least-squares problem.

HyBIT does not accept convergence solely from that estimate. The FGMRES implementation recomputes

r = b - A x
||r||_2

at restart boundaries and before reporting convergence.

For nonzero b, the public report uses

relative_residual = ||b - A x||_2 / ||b||_2

This distinction matters for restarted methods, finite precision, and ill-conditioned or strongly nonnormal problems.

9. Jacobi versus ILU(0) on GeneralSquare

The two current choices make different compromises.

ILU(0) is appropriate for FGMRES because FGMRES does not require the preconditioner to define the SPD inner product required by PCG. That does not mean ILU(0) guarantees good convergence: no-fill LU can be weak, ordering sensitive, or numerically troublesome.

The selective row-relative pivot floor is a factorization safeguard. A nonzero adjusted-pivot count should be treated as diagnostic evidence, not hidden as an ordinary successful factorization.

10. Future SymmetricIndefinite route: MINRES

MatrixProblemClass::SymmetricIndefinite is represented in the API but is not currently routed.

The intended Krylov family is MINRES because it can exploit symmetry without requiring A itself to be positive definite. A future preconditioned MINRES path must preserve the self-adjoint structure required by MINRES; in the standard formulation this normally means using an SPD preconditioner.

A generic nonsymmetric ILU(0) should therefore not simply be reused as the default MINRES preconditioner. Doing so would destroy the very symmetry that makes MINRES attractive.

No current HyBIT call silently falls back from SymmetricIndefinite to PCG or FGMRES.

11. Rectangular systems are a separate problem class

Least-squares systems

min ||A x - b||_2

with rectangular A are not just another GeneralSquare case. Future LSQR/LSMR-class support should have an operator contract that exposes the actions required by the method, including transpose application.

Rectangular systems are currently outside the high-level HyBIT solver.

12. Current compatibility map

This table describes current high-level HyBIT routing, not every mathematically possible combination of Krylov method and preconditioner.

The presence of Minres, Gmres, or Bicgstab variants in reporting enums does not by itself mean that the current high-level solver routes a problem through those algorithms.

13. Failure modes: distinguish three different outcomes

Converged

The solver reached its configured residual criterion. This says nothing by itself about mesh quality, units, boundary conditions, material laws, or other modeling correctness.

MaxIterations

The recurrence remained numerically executable but did not reach tolerance inside the iteration budget. Possible causes include a weak preconditioner, poor restart choice, difficult spectrum/nonnormality, bad scaling, or a problem that does not satisfy the declared model assumptions.

This is not the same as numerical breakdown.

Breakdown

The recurrence encountered a condition that invalidates or prevents the next mathematical step.

Examples in the current implementation include:

• PCG observes non-positive/non-finite r^T M^-1 r;
• PCG observes non-positive/non-finite p^T A p;
• FGMRES encounters a singular/non-finite Hessenberg factor;
• FGMRES Arnoldi can reach a happy breakdown; HyBIT recomputes the true residual and reports convergence only if the actual tolerance is met.

When breakdown occurs, increasing only max_iterations is generally not a meaningful remedy. Recheck the declared matrix class, matrix validity, preconditioner, scaling, and factorization diagnostics.

14. Practical selection workflow

For a new matrix family:

1. Establish the mathematical contract of A.
              |
              v
2. Select the corresponding Krylov family.
              |
              v
3. Start with the simplest compatible preconditioner.
              |
              v
4. Independently verify ||b-Ax|| / ||b||.
              |
              v
5. Inspect iterations, breakdown/status, setup, solve time, and memory.
              |
              v
6. Strengthen one component at a time.
              |
              v
7. For repeated RHS, separate prepared reuse from warm-start effects.

For SPD problems, a PCG breakdown is evidence to investigate the SPD/operator or preconditioner contract, not a reason to disguise the matrix as GeneralSquare.

For GeneralSquare problems, compare Jacobi and ILU(0), inspect ILU pivot adjustments, test ordering sensitivity, then tune restart. Do not infer a universal restart from one matrix family.

For every path, independently recompute the true residual from the original matrix and RHS before treating solver convergence as validated.

HyBIT prepared execution

This document describes analyze -> prepare -> solve-many. HyBIT 0.7.0 provides the SPD and structural prepared paths; develop/0.8.0 additionally provides prepared GeneralSquare FGMRES with Jacobi or opt-in ILU(0).

Analyze

HybitSolver::analyze_csr32 validates the declared problem class, profiles the matrix, resolves backend/execution policy, and records exact structure/value signatures. Analyze does not require an RHS.

Configuration snapshot

Configure the solver before analyze/prepare. Analysis binds the declared problem/execution contract to the matrix signatures. Preparation then copies the relevant solver and preconditioner options into the returned prepared context.

An already-created prepared context is self-contained: later changes to the HybitSolver object do not alter that prepared context. Build a new prepared context when intentionally changing the matrix-dependent solver policy.

Prepare

Preparation creates matrix-dependent reusable state.

Generic SPD

The prepared state may include Jacobi, optional ABTM, reusable PCG workspace, explicit algebraic coarse state, and later learned Hybrid local factors.

GeneralSquare

The prepared state owns exactly one selected preconditioner: Jacobi or ILU(0). It also owns reusable FGMRES workspace. Escalating/budget-aware policies allocate workspace for the maximum permitted restart.

general_square_preconditioner_bytes() reports persistent preconditioner state. For ILU(0), general_square_ilu_adjusted_pivots() reports factor-pivot stabilization count.

Structural SPD

prepare_structural_csr32 additionally owns geometry-dependent aggregation, rigid-body coarse state/factorization, and reusable PCG workspace.

Solve-many

analyze(matrix, policy)
        |
        v
prepare(matrix, analysis)
        |
        +--> solve(rhs_1)
        +--> solve(rhs_2)
        +--> ...

SPD Hybrid can learn local regions/factors on an early difficult RHS and reuse them later.

GeneralSquare builds Jacobi or ILU(0) during prepare and reuses it across all RHS vectors. Its report marks preconditioner_reused = false on solve sequence 1 and true from solve sequence 2 onward, meaning that the same prepared preconditioner is being reused across RHS solves rather than rebuilt.

Report accounting

The first solve in a prepared context charges the stored analysis/prepare time to its SolveReport; later solves report zero for those two setup components.

relative_residual is ||b-Ax||_2 / ||b||_2 for nonzero RHS. This makes solve-many comparisons meaningful even when RHS magnitudes differ.

Exact matrix reuse rule

Prepared contexts validate both CSR structure and exact f64 coefficient bit patterns. Changed values require a new prepare even if sparsity is unchanged.

Future work may separate:

1. same topology + same values -> full reuse;
2. same topology + changed values -> symbolic reuse + numeric refactor;
3. changed topology -> full analyze/prepare.

Only the first case is supported by the high-level prepared context today.

Allocation behavior

Repeated Krylov execution reuses preallocated vectors.

• PCG owns reusable full-size workspace.
• FGMRES owns V, Z, Hessenberg/Givens, and residual/update scratch.
• Local Schwarz owns local scratch.
• Structural preparation owns coarse/geometry state.

First-time diagnostics/factor construction may allocate outside the repeated Krylov inner loop.

Example

let analysis = solver.analyze_csr32(&a)?;
let mut prepared = solver.prepare_csr32(&a, &analysis)?;

let mut x1 = vec![0.0; a.nrows()];
let report1 = prepared.solve(&a, &b1, &mut x1)?;

let mut x2 = vec![0.0; a.nrows()];
let report2 = prepared.solve(&a, &b2, &mut x2)?;

assert_eq!(report1.solve_sequence, 1);
assert_eq!(report2.solve_sequence, 2);
assert!(report2.preconditioner_reused);

When benchmarking reuse, avoid warm-starting from the previous solution unless warm-start behavior is itself the subject of the benchmark.

GeneralSquare FGMRES and ILU(0)

Development documentation for develop/0.8.0, current through F7.

Purpose

MatrixProblemClass::GeneralSquare routes real nonsymmetric square systems through restarted FGMRES. The prepared preconditioner is selected independently from restart policy.

Current choices:

• Jacobi — default;
• ILU(0) — explicit opt-in.

GeneralSquare currently supports ExecutionPolicy::Auto and ExecutionPolicy::Cpu.

For the cross-method rationale—PCG versus FGMRES, preconditioner compatibility, restart complexity, and the future MINRES route—see SOLVER_SELECTION.md.

Why FGMRES

FGMRES stores the Arnoldi basis and the preconditioned basis separately. At Arnoldi step j, conceptually:

v_j                 Arnoldi basis vector
z_j = M_j^-1 v_j    preconditioned basis vector
w   = A z_j

w is orthogonalized against V, and the solution update uses Z:

x = x_0 + Z y

The current high-level path uses a fixed Jacobi or fixed ILU(0), but the separate V/Z representation keeps the recurrence valid for future variable preconditioning.

The implementation uses two-pass modified Gram-Schmidt, Givens rotations, and true residual recomputation at restart boundaries and before accepting convergence.

Restart and memory

Larger restart can reduce restart boundaries or iterations, but increases Arnoldi V/Z storage and orthogonalization work. Larger is therefore not automatically faster.

GeneralSquareOptions controls restart independently from the preconditioner.

Fixed

Default. Default restart is 30.

Escalating

Runs bounded stages and doubles restart up to max_restart. The final maximum restart stage receives the remaining iteration budget.

BudgetAware

Stays inside one FGMRES invocation, observes exact restart-boundary residuals, estimates logarithmic decay over a small window, and grows restart only under sustained projected-budget pressure with a non-improving trend, or stronger emergency pressure.

This is a controller heuristic, not a convergence proof.

Jacobi

GeneralSquare Jacobi permits negative diagonal values but requires a complete, finite, nonzero diagonal.

Advantages are very low setup/storage cost and simple baseline behavior. Its limitation is that it ignores off-diagonal coupling.

ILU(0)

ILU(0) approximates

A ~= L U

without adding fill outside the chosen sparse pattern.

Caller CSR rows need not be sorted or unique. HyBIT privately canonicalizes the preconditioner input:

1. collect each row;
2. sort by column;
3. sum duplicate columns;
4. remove exact-zero off-diagonal entries;
5. verify a diagonal entry.

The original caller matrix remains unchanged and is still used by the sparse operator. Canonicalization belongs only to ILU preparation.

The numeric factorization preserves the canonical sparsity pattern. It performs no row pivoting and introduces no fill.

Factor-pivot stabilization

No-pivot LU/ILU can encounter a zero factor pivot even for a nonsingular matrix. HyBIT records each row scale from the canonical matrix before numeric factor updates, then uses the selective floor

row_scale_i = max_j |a_ij|
floor_i     = 1e-12 * row_scale_i

when a factor pivot magnitude is at or below the floor. Negative small pivots keep negative sign; exact zero receives the positive floor.

The prepared context exposes:

• general_square_ilu_adjusted_pivots();
• general_square_preconditioner_bytes().

A nonzero adjusted-pivot count is diagnostic information. The floor is a factorization safeguard, not a claim that arbitrary near-singular systems are well-conditioned.

Selecting ILU(0)

use hybit::{
    GeneralSquarePreconditionerPolicy, HybitSolver, MatrixProblemClass,
};

let mut solver = HybitSolver::new();
solver.set_problem_class(MatrixProblemClass::GeneralSquare);
solver.set_general_square_preconditioner_policy(
    GeneralSquarePreconditionerPolicy::Ilu0,
);

Selecting ILU(0) does not change restart.

E5 evidence

Synthetic E5 screening compared Jacobi, Ruiz-scaled Jacobi, contiguous dense-LU block Jacobi, and ILU(0). ILU(0) was the strongest tested candidate.

A restart sweep with prepared ILU(0) found restart 3 fastest on all three tested hard families, even though larger restart sometimes used fewer iterations.

The production public-path hard-B cross-check at grid 256 measured:

FGMRES restart          3
iterations              276
true relative residual  9.455610e-9
ILU state               4.238 MiB
FGMRES workspace        4.500 MiB
persistent total        8.738 MiB
adjusted ILU pivots     0

Prepared solves #2 and #3 reproduced the same iteration count and residual while reusing the preconditioner.

These are development regression measurements, not universal performance claims.

F1 ordering evidence

F1 added a benchmark-only Natural-versus-RCM cross-check without changing the production GeneralSquare policy. RCM is formed from the undirected sparsity graph of A + A^T; the benchmark applies the same permutation to rows, columns, and RHS, then maps the computed solution back to the original ordering and recomputes the true residual against the original A and b.

Jacobi is run under both orderings as a control. Because a simultaneous row/column permutation should not materially change diagonal scaling as a mathematical preconditioner, matching Jacobi behavior is used to validate the permutation/verification path before interpreting ILU(0) differences.

Measured fixed-restart-30 results at relative tolerance 1e-8:

For cfd1, ILU(0) solve wall time fell from about 18.0 s to 10.5 s; the measured RCM graph/order/permutation cost was about 64 ms. For raefsky3, solve wall time fell from about 141.6 ms to 51.0 ms, while ordering cost about 41.2 ms. sherman5 demonstrates the opposite regime: Natural ILU(0) already converged in one restart cycle, so RCM added cost without reducing iterations.

All three Natural/RCM ILU(0) comparisons reported zero adjusted pivots. The observed convergence differences therefore cannot be attributed to the row-relative pivot floor.

The corpus also shows that bandwidth reduction is not a sufficient automatic selection signal: sherman5 had the largest relative bandwidth reduction but no ILU(0) iteration reduction. A future automatic policy should therefore consider solve/progress behavior rather than promoting RCM solely from a graph metric.

Goodwin_010 was rejected by the current high-level GeneralSquare contract because at least one structural diagonal entry is missing. F1 deliberately did not insert artificial diagonal values; this matrix is retained as evidence for the later unsuitable-ILU/fallback design.

These measurements use generated b = A * 1 to provide a known solution and an independent forward-error cross-check. They are development evidence for the measured matrices, RHS construction, hardware, and revision rather than a universal ordering rule.

F2 real nonsymmetric corpus evidence

F2 adds preflight classification and a multi-matrix corpus wrapper. Unsupported cases remain recorded evidence instead of aborting the sweep.

Unsupported inputs are not repaired by inserting artificial diagonal values.

For venkat25, generated b=A*1, restart 30, tolerance 1e-8, maximum 5000 iterations produced:

RCM graph/order/permutation cost 52.990 ms. Charging analysis, preparation, ordering, and solve once gives approximately 956.95 ms Natural versus 897.19 ms RCM. Repeated RHS reuse should improve RCM amortization further if the prepared ordering and factors are reused.

Together with F1, the supported real nonsymmetric corpus deliberately contains different ordering regimes: no convergence benefit on sherman5, strong benefit on raefsky3, and moderate benefit on the larger venkat25. Bandwidth reduction is therefore telemetry, not an automatic selection rule.

F3 prepared multi-RHS ordering evidence

F3 measures whether the F1/F2 ordering behavior survives prepared solve-many reuse with different right-hand sides. Natural and RCM ILU(0) states are each prepared once and reused for five deterministic RHS vectors (ones, alternating, then deterministic hashed vectors). Every solve starts from a zero initial guess. The RCM path includes per-RHS P b construction and solution unpermutation in the end-to-end accounting.

The benchmark verifies the final residual in the original system after both timed solve paths complete. From solve sequence 2 onward it also requires the prepared context to report preconditioner reuse.

Five repeated 5-RHS runs produced:

For raefsky3, the same ILU(0) ordering advantage seen with b=A*1 remains strong across all five RHS families. The deterministic iteration totals are 225 Natural versus 69 RCM over the five RHS vectors. The median cumulative permutation/unpermutation cost is only about 0.257 ms, so end-to-end behavior is dominated by the Krylov/triangular-apply savings.

venkat25 behaves differently. Its deterministic five-RHS iteration total is 737 Natural versus 692 RCM, but most of that difference comes from the first two RHS vectors. Later hashed RHS vectors differ by only a few iterations and individual RCM wall times can match or exceed Natural. Across five repeated process runs the median end-to-end advantage is only about 1.57%, despite RCM being lower in all five paired totals. This is treated as near-parity, timing-sensitive evidence rather than a robust performance win.

The result weakens two potential automatic-selection signals. Structural bandwidth alone was already insufficient after F1/F2; F3 additionally shows that a strong result on one RHS does not necessarily predict the ordering benefit for later RHS vectors on the same matrix. A future ordering policy should therefore use progress/solve evidence conservatively rather than promoting RCM from bandwidth or one isolated RHS result.

F4 serial triangular-apply and dependency evidence

F4 separates ILU(0) triangular application cost from the surrounding FGMRES solve and compares that kernel with serial CSR SpMV. The retained benchmark also reconstructs the canonical ILU(0) pattern and profiles forward/backward dependency levels without changing the production solver.

The measured serial kernel medians were:

The three matrices therefore expose different mechanisms. sherman5 gets a faster triangular kernel from RCM but no measured F1 iteration reduction. raefsky3 gets essentially no triangular-kernel speedup, so its strong RCM benefit is dominated by improved ILU(0) effectiveness. venkat25 gets a faster SpMV but a slower triangular application, partially offsetting its modest RCM iteration reduction.

The dependency profile also shows why bandwidth alone is insufficient:

• sherman5 RCM changes forward/backward levels from 39/66 to 52/70 while reducing dependency-distance p95 from 1105 to 117.
• raefsky3 RCM changes 1176/1176 levels to 1592/1592 and reduces structural average parallelism from about 18.03 to 13.32.
• venkat25 RCM changes 4176/4176 levels to 1700/1700 and raises structural average parallelism from about 14.95 to 36.72, but dependency-distance median grows from 117 to about 711/712 and p95 from 1948 to 2348.

A follow-up experimental level-scheduled Rayon prototype produced numerically identical results, but one parallel dispatch/barrier per triangular level was too fine grained. A width-threshold hybrid did not recover the overhead: every measured case that actually invoked Rayon was slower than the canonical serial row-order apply. That prototype is intentionally not retained in the production code.

The current conclusion is to keep canonical ILU(0) triangular application serial. Future parallel work, if justified by a larger workload, should use a coarser persistent-worker, task-DAG, or superlevel design rather than the rejected per-level Rayon schedule.

Limitations

• ILU(0) is ordering-sensitive.
• Current L/U triangular application is serial.
• There is no ILUT, drop tolerance, level-of-fill, threshold pivoting, or internal reordering.
• F2 validates three supported real nonsymmetric ordering cases up to 62424 unknowns, but the corpus is not exhaustive; five additional screened matrices expose the current complete-diagonal contract as a major applicability boundary.
• GeneralSquare is not yet routed through resident Rayon or GPU execution.
• GeneralSquare selection is not yet exposed through the C ABI configuration surface.

Recommended evaluation

For a new GeneralSquare family:

1. run Jacobi baseline;
2. verify true residual;
3. opt into ILU(0);
4. record factor bytes and adjusted pivot count;
5. sweep a small set of restart values;
6. compare natural and at least one meaningful reordering;
7. benchmark repeated RHS reuse;
8. validate representative real matrices before changing automatic policy.

F5 unsuitable-ILU fallback

F5 separates solver applicability from ILU(0) applicability.

The existing GeneralSquarePreconditionerPolicy::Ilu0 remains strict. A structurally missing diagonal still rejects explicit strict ILU(0). The new explicit GeneralSquarePreconditionerPolicy::Ilu0Fallback first attempts the same canonical ILU(0) preparation and falls back to a prepared Identity preconditioner only for HybitError::MissingDiagonal. Other ILU preparation errors remain errors and are not silently hidden.

The prepared context exposes the effective preconditioner kind and whether the structural fallback was used. A solve using the Identity fallback reports PreconditionerKind::None.

The real missing-diagonal corpus from F2 was rechecked with this policy. All five matrices (Goodwin_010, Goodwin_023, Goodwin_030, goodwin, and rma10) retained strict-ILU rejection while Ilu0Fallback reached prepared state with fallback_used=true.

A retained bounded-solve harness then tested the Identity fallback for 300 FGMRES iterations at restart 30 and tolerance 1e-8 using deterministic manufactured right-hand sides. The verified relative residuals were:

Identity is therefore retained as a correctness/safety fallback, not promoted as a claim of strong numerical preconditioning.

F5 also explored several benchmark-only stronger fallbacks and rejected them from the retained implementation:

• inserting row-relative synthetic diagonal values before ILU(0) helped the three Goodwin_0xx matrices but degraded goodwin and especially rma10;
• arbitrary structural perfect matching existed for all five matrices but the resulting column-permuted ILU(0) was numerically unstable;
• row-relative bottleneck matching removed the weakest matched entries but all five cases still broke down;
• alternating max-norm row/column equilibration drove row and column norms near one after eight sweeps, yet the matching-ILU triangular solve still broke down. goodwin and rma10 did so with zero adjusted pivots.

These negative experiments delimit the current canonical ILU(0) operating envelope: missing-diagonal systems need more than structural completion, matching, or simple equilibration for a robust stronger preconditioner. The rejected F5c-F5f experiment programs are not retained in the repository.

The next GeneralSquare work is policy selection from residual/progress evidence, not another ad-hoc missing-diagonal transformation. Drop-tolerance/fill or pivoting methods may be considered later as separate preconditioners now that the ILU(0) envelope is documented.

F6 ordering-selection evidence

F6 tested whether Natural -> RCM ILU(0) ordering can be selected from cheap runtime evidence without promoting RCM unconditionally.

F6a: Natural-only progress is insufficient

F6a collected exact FGMRES restart-boundary residuals under Natural ILU(0) and compared them with the eventual Natural-versus-RCM result.

Natural-only early progress did not separate the measured regimes. For example, sherman5 and raefsky3 both showed strong first-cycle logarithmic residual decay, but RCM provided essentially no convergence benefit on sherman5 and a large benefit on raefsky3. venkat25 showed much slower Natural early decay yet only modest, RHS-dependent RCM benefit.

A first-cycle projected-total-iteration extrapolation was also unreliable because convergence rates changed after restart boundaries.

F6b: paired short probes are more discriminating

F6b prepared both Natural and RCM ILU(0) states and compared equal-length 4-, 8-, and 16-iteration FGMRES probes from x=0.

The key benchmark signal is

paired_ratio =
    (RCM residual / initial residual)
    /
    (Natural residual / initial residual)

after the same number of probe iterations.

The 4-iteration ratio separated several important measured cases:

The additional capped cases are classification evidence only; they are not claimed as converged performance wins.

cfd1 also showed why the paired signal must not be interpreted as a per-RHS winner predictor. RCM was much faster for RHS 1, 2, and 4 but slower for RHS 3 and 5. The first-RHS paired probe is therefore best interpreted as a candidate signal for a prepared solve-many workload, not as a guarantee for every subsequent right-hand side.

F6c: selector cost must be charged

F6c replayed the candidate decision while charging Natural preparation, RCM ordering and preparation, both short probes, and the selected full solves.

With a four-iteration paired probe and a provisional threshold of 1.0:

• cfd1 selected RCM with paired ratio 0.945698. Across five repeats, the median policy/Natural cost ratios were 0.591, 0.606, 0.733, 0.731, and 0.798 through RHS counts 1 through 5.
• thermal1 selected RCM with paired ratio 0.579999. The corresponding median ratios were 0.780, 0.625, 0.616, 0.622, and 0.609.
• raefsky3 had already shown paired ratio 0.249599; with the earlier replay it became clearly profitable by the third RHS and reached about 0.668 through five RHS vectors.
• venkat25 correctly retained Natural with paired ratio 2.170262, but paying for the rejected alternate state still left about 3.8% overhead after five RHS vectors.
• a provisional small-row guard avoided this overhead on sherman5, but row count is a benchmark heuristic rather than a validated production criterion.

F6 conclusion

The evidence does not justify an unconditional production Natural -> RCM promotion.

A useful experimental solve-many selector now has two evidence-backed pieces:

1. skip RCM when it does not reduce structural bandwidth;
2. when RCM does reduce bandwidth, a four-iteration paired residual ratio below one is a promising numerical signal for an RCM candidate.

However, the paired probe requires constructing and probing the alternate RCM state even when Natural is ultimately retained. venkat25 demonstrates a real selection-overhead regression, and cfd1 demonstrates right-hand-side sensitivity after the matrix-level decision.

Therefore F6 retains the progress, paired-probe, and policy-replay harnesses as development evidence but does not change the production ordering/default policy. A later explicit solve-many ordering policy may reuse these signals if the caller can provide an expected reuse horizon or another cost-aware trigger.

F7 Jacobi -> ILU(0) promotion evidence

F7 tested whether GeneralSquare can safely promote from the default Jacobi preconditioner to canonical Natural-order ILU(0) automatically.

The answer for the current 0.8 development evidence is no: ILU(0) is extremely valuable on several matrices, but no cheap signal tested in F7 separated those wins from the measured regressions without important false positives or false negatives.

F7a: complete prepared Jacobi-versus-ILU(0) solves

F7a charged analysis, preparation, and solve wall time for both prepared preconditioners on deterministic right-hand sides.

Representative RHS1 results at restart 30 and tolerance 1e-8:

The same canonical ILU(0) implementation therefore ranges from a decisive robustness/performance improvement to a large setup/apply penalty with worse convergence.

Simple row density cannot classify these cases. For example raefsky3 has roughly 70 nonzeros per row and strongly benefits from ILU(0), whereas cant and x104 are of comparable row density and regress.

F7b: one-apply approximate-inverse probes

F7b measured static structural/cost signals and a deterministic one-apply quality probe

defect(M) = ||M^-1 A x - x|| / ||x||

for Jacobi and ILU(0).

This signal is useful as a catastrophic-factor-quality diagnostic but is not a positive promotion selector:

• raefsky3 and venkat25 strongly benefit from ILU(0) in complete FGMRES despite median ILU/Jacobi defect ratios greater than one.
• boneS01 has a favorable median defect ratio below one yet remains slower with ILU(0) and fails to reach tolerance in the complete run.
• nd3k, cant, s3dkq4m2, and x104 show strong negative signals, including very large defect ratios or ILU work/setup cost.

Accordingly F7b can veto obviously pathological factors in experiments, but no production threshold is assigned from the current corpus.

F7c: short paired Krylov probes

F7c compared actual Jacobi-FGMRES and ILU(0)-FGMRES progress from the same x=0 state after 4, 8, and 16 iterations.

The paired Krylov signal is more informative than the one-apply defect, but it still is not sufficient for promotion:

• four iterations are too short: raefsky3 initially gives an ILU/Jacobi residual ratio above one even though ILU(0) is a major complete-solve win;
• venkat25, cfd1, and thermal1 show clear early ILU residual advantages;
• nd3k, s3dkq4m2, and x104 are rejected clearly;
• cant and boneS01 are false positives: ILU(0) has a smaller early residual while the complete wall time is substantially worse.

F7d: restart-boundary replay

F7d replayed the complete prepared comparison with iteration budgets of 30, 60, 120, and 240, corresponding to one, two, four, and eight restart-30 cycles.

Restart-boundary telemetry explains several F7c false positives but does not produce a universal selector.

cant is the clearest reversal. At 30 iterations ILU(0) has the smaller residual, but by 60 iterations Jacobi is already better; the ILU residual then remains worse through 120 and 240 iterations.

boneS01 is the opposite difficulty. ILU(0) keeps the smaller residual through 240 iterations, but its higher setup and per-iteration cost keeps total wall time substantially above Jacobi. At 240 iterations the ILU/Jacobi total-time ratio is still about 1.81.

raefsky3 and venkat25 show why a simple early cost-normalized rule is also unsafe. Their ILU setup/apply cost can make short truncated comparisons look unfavorable even though ILU(0) later reaches the requested tolerance far before Jacobi, or succeeds when Jacobi reaches the 5000-iteration cap.

The convergence rate is therefore not stationary enough across restarts to project time-to-tolerance safely from a small fixed observation window.

F7 conclusion

HyBIT 0.8 keeps GeneralSquare Jacobi as the default and canonical ILU(0) as an explicit opt-in policy.

F7 does not add automatic Jacobi -> ILU(0) promotion.

This is a deliberate robustness decision rather than a statement that Jacobi is generally better. ILU(0) is dramatically superior on part of the measured corpus and can turn a 5000-iteration Jacobi failure into convergence. The problem is safe automatic classification: the tested static, factor-quality, short-Krylov, and restart-boundary signals all leave meaningful counterexamples.

The retained F7a/F7b/F7c harnesses are development evidence for future cost-aware or application-informed policy work. A future automatic policy should require broader held-out validation and may need caller-provided solve horizon, memory budget, or an explicit willingness to pay for exploratory factorization/probing.

Structural Auto API (0.8)

HyBIT’s generic solve_csr32 path remains geometry-free and backward compatible. For 3-D structural SPD systems, callers can provide node coordinates explicitly:

use hybit::{
    HybitSolver, RigidBodyAggregation, SolverOptions, StructuralOptions,
    StructuralPcgVectorPolicy, StructuralPreconditionerPolicy, StructuralSpmvPolicy,
};

let mut solver = HybitSolver::new();
solver.set_options(SolverOptions {
    relative_tolerance: 1.0e-8,
    absolute_tolerance: 0.0,
    max_iterations: 3000,
})?;
solver.set_structural_options(StructuralOptions {
    target_coarse_dimension: 1536,
    aggregation: RigidBodyAggregation::Auto,
    spmv_policy: StructuralSpmvPolicy::Auto,
    preconditioner_policy: StructuralPreconditionerPolicy::Auto,
    pcg_vector_policy: StructuralPcgVectorPolicy::Auto,
})?;

let report = solver.solve_structural_csr32(&a, &coordinates, &b, &mut x)?;

For solve-many workloads, prepare once and reuse the rigid-body coarse factor:

let analysis = solver.analyze_csr32(&a)?;
let mut prepared = solver.prepare_structural_csr32(&a, &analysis, &coordinates)?;
let first = prepared.solve(&a, &b1, &mut x1)?;
let second = prepared.solve(&a, &b2, &mut x2)?;

HybitPreparedStructuralSystem exposes the selected aggregate size, aggregate count, coarse dimension, factor storage, geometry storage, and workspace size.

The structural path currently assumes exactly three displacement DOFs per node and uses six rigid-body modes per aggregate. RigidBodyAggregation::Auto is the default: it tries deterministic graph-connected aggregation first and falls back to the contiguous RCM-order baseline if the graph coarse Cholesky is numerically singular or if the structural node graph contains disconnected components too small to support a six-mode rigid-body aggregate. RigidBodyAggregation::Contiguous preserves the fixed-width baseline, while RigidBodyAggregation::Graph requests graph aggregation without fallback. Coordinate ordering must match the CSR DOF ordering exactly: node i maps to rows 3*i, 3*i+1, 3*i+2.

Measuring prepared reuse

bench-fem-structural-prepared.ps1 prepares the structural context once and then solves the same RHS repeatedly from a fresh zero initial guess. This deliberately prevents warm-start effects from being confused with preconditioner reuse. The first report charges analysis and preparation; later reports must show preconditioner_reused = true with zero analysis/prepare charge.

Example:

.\benchmarks\scripts\bench-fem-structural-prepared.ps1 K.mtx `
  -Coordinates K.coords `
  -Rhs b.txt `
  -TargetCoarseDimension 1536 `
  -Tolerance 1e-8 `
  -MaxIterations 3000 `
  -Repeats 2

Graph aggregation experiment (0.6-r10)

The graph mode keeps the same target coarse dimension and six rigid-body modes, but forms connected node regions from 3x3 block sparsity instead of slicing the RCM node order into fixed-width chunks. This is intended as an A/B experiment before changing the default Structural Auto policy.

.\benchmarks\scripts\bench-fem-structural-auto.ps1 K.mtx `
  -Coordinates K.coords `
  -Rhs b.txt `
  -TargetCoarseDimension 1536 `
  -Aggregation graph `
  -Tolerance 1e-8 `
  -MaxIterations 3000

The benchmark reports actual aggregate count and minimum/maximum aggregate sizes so changes in coarse dimension are visible rather than hidden.

Graph aggregation robustness (0.6-r11)

The initial capped-BFS graph experiment could leave tiny remainder islands after nearby nodes had already been assigned. A six-mode rigid-body basis can become rank deficient or severely ill-conditioned on such islands. r11 therefore merges graph aggregates smaller than one quarter of the requested aggregate size (never below three nodes) into the neighboring aggregate with the strongest edge connection before forming Z^T A Z. No diagonal regularization is added: failure of the coarse Cholesky still indicates a genuine basis/geometry problem rather than being hidden by a shift.

Structural Auto graph policy (0.6-r12)

The L-angle physical-load development benchmark showed that graph-connected aggregation was substantially more effective than contiguous RCM chunks at essentially the same coarse dimension: 220 versus 1370 PCG iterations at a relative tolerance of 1e-8, with coarse dimensions 1398 and 1404 respectively. This is a single development benchmark, not a general performance guarantee.

For that reason r12 makes RigidBodyAggregation::Auto the Structural Auto default. Auto attempts the graph coarse space first. If coarse Cholesky reports a numerical breakdown, or if graph topology leaves a disconnected component with fewer than three nodes, Auto retries with contiguous aggregation. Explicit Graph remains strict so A/B benchmarks still expose graph-basis failures rather than silently hiding them.

Packed coarse Cholesky storage (0.6-r17)

The structural rigid-body path now stores the dense Cholesky lower factor in row-wise packed triangular form. The Galerkin coarse operator, aggregation, rigid-body basis, pivot test, and PCG algorithm are unchanged; only factor storage and triangular-solve addressing change. A coarse dimension m uses m(m+1)/2 stored f64 factor values instead of m^2, with row-start offsets precomputed once during preparation.

This is intentionally an implementation optimization rather than a policy change. Re-run the Graph target-coarse-dimension points at 1536 and 3072 to measure whether the reduced factor footprint improves wall time as well as memory use on the target CPU.

Fine-grid SpMV policy (0.6-r20)

Structural PCG can now select the validated Rayon-parallel CSR operator through StructuralSpmvPolicy. Auto keeps small matrices serial and switches CSR32 systems with at least 1,000,000 nonzeros to parallel row execution. Serial and Parallel remain available for controlled benchmarking; explicit Parallel requires the CSR32 backend. The parallel implementation preserves the same CSR storage and arithmetic ordering within each row, so iteration counts should remain unchanged.

The Rayon pool size is intentionally external to the solver policy. Rust callers can continue to use RAYON_NUM_THREADS. C/C++/Fortran hosts additionally have the 0.8 thread-interoperability API: explicit hybit_set_num_threads(n) takes precedence, followed by RAYON_NUM_THREADS, then OMP_NUM_THREADS, then the Rayon default. OpenMP API users should synchronize omp_get_max_threads() into HyBIT before creating the first solver. See Threading and host-runtime interoperability. On the development L-angle case on a Ryzen 7 7800X3D, 4 threads was the best tested point: CSR SpMV improved from about 13.24 ms to 5.09 ms and the 220-iteration structural PCG solve from about 4.44 s to 2.73 s. This machine-specific result is not used as a hard-coded library thread count.

0.6 structural execution policy (development)

Structural Auto independently resolves StructuralSpmvPolicy, StructuralPreconditionerPolicy, and StructuralPcgVectorPolicy. Large CSR structural systems can therefore use parallel CSR SpMV, parallel 3x3 block-Jacobi/restriction/prolongation, and parallel/fused PCG vector kernels together. Small systems remain serial to avoid Rayon overhead. The packed dense coarse triangular solve remains serial. Rayon thread count is controlled externally (for example RAYON_NUM_THREADS).

PCG dense-vector policy (0.6-r25)

Structural solves expose StructuralPcgVectorPolicy::{Auto, Serial, Parallel} independently from the sparse SpMV and rigid-body preconditioner policies. The parallel path uses the shared Rayon pool for dot/norm reductions, search-direction updates, and a fused x += alpha*p, r -= alpha*Ap, ||r|| traversal.

Auto is deliberately conservative: it selects the parallel/fused path only for at least 131072 unknowns and at least 4 Rayon workers. Explicit Serial remains available for low-thread-count execution and controlled A/B tests.

At the r25 production checkpoint, the physical-load 358065-DOF / 28239653-nnz L-angle case at 8 Rayon workers resolved to Graph aggregation, Parallel SpMV, Parallel preconditioning, and Parallel PCG vectors. It used 233 aggregates and a 1398-dimensional coarse space, converged in 220 iterations, and independently verified a relative residual of 9.378557e-9. Measured solve time was 1.726 s; analysis+prepare+solve was 2.797 s on the Ryzen 7 7800X3D development machine. These timings are a reproducible development reference, not a cross-machine performance guarantee.

Threading and host-runtime interoperability

HyBIT uses Rayon for validated CPU-parallel execution. C, C++, and Fortran applications frequently already use OpenMP, and FEM hosts such as FrontISTR may combine MPI ranks with OpenMP threads. HyBIT 0.8 therefore treats thread-count selection as an explicit process-level interoperability contract.

Thread-count precedence

Before the first HyBIT solver is created, the effective Rayon worker count is selected in this order:

1. hybit_set_num_threads(n)
2. RAYON_NUM_THREADS
3. OMP_NUM_THREADS
4. Rayon default

For an OpenMP nesting list such as:

OMP_NUM_THREADS=8,4,2

HyBIT uses the first value (8) for its Rayon worker pool.

RAYON_NUM_THREADS intentionally takes precedence over OMP_NUM_THREADS. This permits a host application to keep its OpenMP team size and HyBIT’s Rayon pool size different when that is deliberate.

OpenMP API settings

An environment-variable fallback cannot observe a later call such as:

omp_set_num_threads(8);

For applications that select the OpenMP team size through the OpenMP API, call HyBIT’s thread API before creating the first solver:

omp_set_num_threads(8);
int rc = hybit_set_num_threads((uint32_t)omp_get_max_threads());

When the C/C++ translation unit is compiled with OpenMP, hybit.h also exposes the header-only convenience bridge:

hybit_sync_openmp_threads();

and the C++ wrapper exposes:

hybit::sync_openmp_threads();

Fortran can use the normal OpenMP module and the ISO_C_BINDING wrapper:

use omp_lib
use hybit
use, intrinsic :: iso_c_binding

integer(c_int) :: rc

call omp_set_num_threads(8)
rc = hybit_set_num_threads(int(omp_get_max_threads(), c_int32_t))

hybit_num_threads() reports the effective Rayon worker count and is intended for diagnostics. Because querying it may initialize Rayon’s global pool, set the desired policy first.

Process-wide lifetime

Rayon’s global worker pool is process-wide and cannot be resized after initialization. HyBIT therefore does not pretend that thread count is a solver-local property.

Repeated calls to hybit_set_num_threads(n) with the same value are accepted. A request for a different value after initialization is rejected.

Configure the thread count before:

• hybit_solver_create;
• any call to hybit_num_threads;
• any other Rust component in the process that may initialize Rayon’s global pool.

If another Rust library initializes Rayon first with a different size, HyBIT cannot resize that pool.

OpenMP runtime linkage

HyBIT does not link directly to libgomp, LLVM libomp, Intel iomp, or another OpenMP runtime merely to discover a thread count. This avoids silently introducing a second OpenMP runtime into a host process.

The OpenMP helper in hybit.h is header-only and calls the OpenMP API selected by the host application’s own compiler/linker.

Avoid nested OpenMP x Rayon oversubscription

Do not invoke the same parallel HyBIT solve simultaneously from every worker of an active OpenMP team.

A typical integration should call one HyBIT solve from one controlling host thread, outside an active OpenMP parallel region when practical. Application threads should be synchronized before and after the call as required by the host program.

Otherwise an OpenMP team of T threads can each enter a Rayon computation with up to T workers, causing severe oversubscription.

Prepared HyBIT contexts are not intended to be concurrently mutated by multiple host threads. Internal Rayon parallelism is separate from concurrent host calls.

MPI + OpenMP hosts

HyBIT 0.8 does not yet implement a distributed MPI solver, but it can be linked into an MPI/OpenMP application.

For a host with R MPI ranks per node and T CPU threads allocated to each rank, the normal integration is:

MPI rank 0 -> one HyBIT call -> Rayon pool of T workers
MPI rank 1 -> one HyBIT call -> Rayon pool of T workers
...
MPI rank R-1 -> one HyBIT call -> Rayon pool of T workers

Size the job allocation so the node has enough CPUs for approximately R * T software workers. The thread count is per process/rank.

If the MPI launcher sets a CPU set or process affinity for each rank, Rayon workers normally execute within the process scheduling constraints established by the OS/launcher. HyBIT does not translate OpenMP OMP_PLACES or OMP_PROC_BIND policy into Rayon-specific affinity rules. Thread-count synchronization and thread affinity are separate concerns.

FrontISTR-style integration

For a FrontISTR-like MPI/OpenMP host, two straightforward configurations are:

Environment-controlled:

OMP_NUM_THREADS=8
RAYON_NUM_THREADS unset

HyBIT then uses 8 Rayon workers when its pool is first initialized.

API-controlled:

omp_set_num_threads(8);
hybit_set_num_threads((uint32_t)omp_get_max_threads());

The second form is preferred when the application changes the OpenMP team size through API calls rather than through the environment.

If RAYON_NUM_THREADS is set, it overrides the OMP_NUM_THREADS fallback.

Performance policy

A matching thread count does not imply that every sparse workload should execute in parallel. HyBIT’s validated execution policies can retain serial kernels for small problems where Rayon scheduling overhead dominates.

G7 also showed that one universal hardware-independent nnz threshold is not sufficient for all prepared local operators. Explicit serial, full-Rayon, and task-limited prepared execution therefore remain available without embedding a machine-specific automatic selector.

Thread count controls the size of the available worker pool; it does not force every operation to use every worker.

SPD Hybrid preconditioner — mathematical notes

This document covers the Hybrid preconditioner used by HyBIT’s MatrixProblemClass::Spd / PCG route. It does not describe the GeneralSquare FGMRES path.

For the matrix-class and Krylov-method selection rules, see SOLVER_SELECTION.md.

Let the global SPD matrix be A. The controller starts with the configured base SPD preconditioner: algebraic two-level coarse when explicitly enabled, otherwise Jacobi. The initial probe_iterations segment is only a controller boundary; if the preconditioner is unchanged, HyBIT resumes the same PCG recurrence. When that stage shows poor progress, the policy identifies hard-core DOF sets and grows each set independently through the ABTM graph by a configurable number of halo layers. For local region H_k, HyBIT forms the principal matrix

A_k = A[H_k, H_k].

Because a principal submatrix of an SPD matrix is SPD, exact Cholesky is mathematically valid in exact arithmetic. The implementation additionally checks numerical symmetry and positive pivots.

Symmetric overlap weighting

Regions may overlap. Let m_i be the number of local regions containing global DOF i. Within every local region containing i, HyBIT uses

w_i = 1 / sqrt(m_i).

For restriction operator R_k and diagonal local weight W_k, the local correction is

R_k^T W_k A_k^-1 W_k R_k.

Each such term is symmetric positive semidefinite. Jacobi is retained only for global DOFs with zero local multiplicity. Consequently every global DOF is covered by at least one positive contribution while a single fully covering local block reduces to the exact local inverse without double-counting Jacobi.

The conceptual preconditioner is

M^-1 = J_uncovered + sum_k R_k^T W_k A_k^-1 W_k R_k.

This construction is intended to remain compatible with standard PCG assumptions. HyBIT restarts PCG only when selective-direct strengthening changes the preconditioner. Controller/telemetry boundaries with an unchanged preconditioner retain the same resumable PCG session and therefore preserve conjugacy.

Bounded memory

Dense local Cholesky factors scale quadratically in local region order. HyBIT therefore caps each expanded local factor and caps the number of local regions. The current Hybrid defaults are 128 DOFs per local region and at most 8 local regions.

The report exposes both total local factor DOFs (overlap counted per factor) and unique covered DOFs, plus an estimate of bytes actually owned by local factors, indices, overlap weights and multiplicity metadata.

Limitations

The current hard-region detector is heuristic. It combines matrix coupling/scale diagnostics with residual seeds; it is not an estimator of the global condition number. Later releases should incorporate KrylovScope/Lanczos/Ritz information and reusable factor policies.

HyBIT architecture

HyBIT 0.6 retains the staged-PCG safety and reusable prepared execution model from 0.5, and adds an explicit geometry-aware structural path. The generic solve_csr32 behavior remains separate from structural policy selection.

0.8 execution-layer direction

The 0.8 development line separates three concepts that were previously easy to conflate:

• MatrixBackend: sparse/operator representation such as CSR32, ABTM, or matrix-free;
• ExecutionTarget: where Krylov work runs (Cpu today, Gpu as the resident-device target);
• MatrixProblemClass: the mathematical contract (Spd, SymmetricIndefinite, or GeneralSquare).

MatrixProblemClass is declarative. HyBIT does not infer positive definiteness from a positive diagonal or other weak structural checks.

hybit-krylov now contains an experimental KrylovExecutionBackend boundary plus ResidentPcgWorkspace. The established pcg_with_workspace implementation remains the production CPU path in 0.8-a1. A CPU reference execution backend runs the same PCG recurrence through resident vectors and is regression-tested against the established path.

The resident workspace owns solution, RHS, and Krylov scratch vectors. A future CUDA/CubeCL backend can therefore upload matrix/preconditioner state during prepare, upload one RHS and initial solution per solve, keep Krylov vectors resident for the complete iteration loop, and download only the final solution/report data. This avoids designing a GPU backend around per-kernel host/device transfers.

The same low-level execution boundary is intentionally broader than PCG so later MINRES and flexible GMRES work can reuse device vector operations. The automatic solver/controller remains real SPD + PCG in this checkpoint; no existing call is silently routed to GPU or to a non-PCG method.

0.8-a2 control plane

HybitSolver now exposes independent ExecutionPolicy and MatrixProblemClass controls. ExecutionPolicy::Auto and ExecutionPolicy::Cpu resolve to the established CPU implementation. ExecutionPolicy::Gpu is recognized but rejected before preparation until a resident GPU backend is implemented and validated.

Likewise, MatrixProblemClass::Spd remains the only executable automatic-solver class in this checkpoint. SymmetricIndefinite and GeneralSquare are represented explicitly so future MINRES and FGMRES/BiCGStab work does not require another API redesign, but neither class is silently routed through PCG.

Analyze records the resolved execution target and declared problem class, and prepare verifies that those contracts have not changed. This keeps prepared-state reuse deterministic as additional execution targets and Krylov methods are introduced.

0.8-a3 resident CPU cross-check

ExecutionPolicy::CpuResident connects the generic HybitSolver prepare/solve path to CpuKrylovExecution and ResidentPcgWorkspace for fixed-Jacobi PCG. This is an explicit validation policy rather than the default: Auto and Cpu continue to use the established staged PCG implementation.

The a3 resident path deliberately requires HybridOptions.enabled = false. This keeps the first production-level cross-check focused on the fixed-preconditioner recurrence that will also form the first CUDA PoC. Public tests compare legacy and resident CPU solves for convergence status, iteration count, residual, and solution, and verify solve-many reuse of the resident vectors.

For this validation checkpoint the prepared object retains both the established five-vector PCG workspace and the seven-vector resident workspace, so the reported Krylov workspace bytes include both. Once resident execution becomes the normal prepared representation, the duplicate legacy workspace can be removed.

0.8-b1 Rayon resident backend

RayonKrylovExecution adds a second CPU implementation of KrylovExecutionBackend without changing the PCG recurrence. It reuses the already-validated parallel/fused dense-vector kernels from the structural CPU path for dot/norm reductions, fused solution/residual updates, and search-direction updates. Generic axpy/scale operations use the same chunked shared Rayon pool.

Sparse operator and preconditioner application remain delegated through the existing traits. Rayon vector execution can therefore compose independently with serial or parallel CSR and preconditioner implementations; this checkpoint does not change HybitSolver, ExecutionPolicy::CpuResident, or production defaults.

The resident Rayon recommendation initially mirrors the validated Structural Auto vector crossover: at least 131072 unknowns and at least four shared Rayon workers. This is a CPU-vector scheduling heuristic, not a new ExecutionTarget.

GPU work is intentionally deferred until CubeCL 0.11 reaches a stable release. MPI remains a later distributed-memory layer after single-node CPU parallel execution is consolidated.

0.8-b2 resident Rayon CSR execution

ExecutionPolicy::CpuResidentRayon is an explicit validation policy that combines RayonKrylovExecution with ParallelCsr32Operator while retaining the existing serial Jacobi preconditioner. It therefore isolates the two already-validated Rayon layers: CSR row-parallel SpMV and resident dense-vector kernels.

CpuResident remains the serial resident comparison path, while Auto and Cpu remain unchanged. CpuResidentRayon requires HybridOptions.enabled = false and the CSR32 matrix backend; forced ABTM is rejected rather than silently dropping parallel CSR execution.

The resident workspace is still reused across right-hand sides. This checkpoint does not yet parallelize generic Jacobi application and does not route any default policy into the resident Rayon path.

0.8-b3 resident Rayon Jacobi execution

ParallelJacobiPreconditioner is a read-only Rayon view over the already prepared inverse diagonal. It adds no persistent numerical storage and leaves the ordinary JacobiPreconditioner trait implementation serial for low-overhead and A/B use.

ExecutionPolicy::CpuResidentRayonJacobi composes three explicit CPU layers: ParallelCsr32Operator, RayonKrylovExecution, and ParallelJacobiPreconditioner. The b2 CpuResidentRayon policy remains available with serial Jacobi, so serial-all, parallel-CSR+vectors, and parallel-CSR+vectors+Jacobi can be compared in one build.

No automatic policy is changed in this checkpoint.

0.8-b4 resident structural cross-check

HybitPreparedStructuralSystem::solve_resident_rayon adds an explicit validation path for the fully parallel structural configuration. It reuses the already prepared RigidBodyTwoLevelBlockJacobiPreconditioner, ParallelCsr32Operator, and ParallelRigidBodyTwoLevelPreconditioner, but runs PCG through RayonKrylovExecution.

The seven-vector resident workspace is allocated lazily on the first resident structural solve and reused for later right-hand sides. Ordinary solve remains the production path and pays no resident-workspace memory cost unless the validation method is called.

This checkpoint intentionally requires effective Parallel policies for SpMV, structural preconditioning, and PCG vectors. It is a same-preconditioner recurrence cross-check, not a new automatic policy.

0.8-e1 flexible GMRES kernel

hybit-krylov now contains a standalone restarted FGMRES implementation for real square operators. It is intentionally not wired into HybitSolver yet, so SPD/PCG and Structural Auto behavior remain unchanged.

FlexiblePreconditioner receives the global Krylov iteration and a mutable receiver. Existing fixed Preconditioner implementations automatically satisfy this interface, while a future adaptive HyBIT preconditioner can change its approximation between Arnoldi steps without invalidating the FGMRES recurrence.

FgmresWorkspace stores the Arnoldi basis V, independently preconditioned basis Z, Hessenberg/Givens state, and full-size scratch vectors for solve-many reuse. The first implementation uses two-pass modified Gram-Schmidt and recomputes the true residual at every restart and before accepting convergence.

This checkpoint deliberately stops below the automatic solver dispatch layer. MatrixProblemClass::GeneralSquare remains rejected by HybitSolver until a later checkpoint explicitly connects the validated FGMRES kernel and defines the nonsymmetric preconditioner policy.

0.8-e2 prepared GeneralSquare FGMRES

An explicitly declared MatrixProblemClass::GeneralSquare now routes the generic prepared solver through restarted FGMRES with fixed diagonal Jacobi. GeneralSquareOptions::restart controls the Arnoldi restart dimension and defaults to 30. Solve reports identify this path explicitly as SolverKind::Fgmres rather than collapsing it into the ordinary GMRES family.

The E2 checkpoint accepts negative diagonal entries for the general-square Jacobi scaling, while still requiring a complete finite nonzero diagonal. The established SPD Jacobi constructor keeps its positive-diagonal requirement unchanged.

GeneralSquare currently supports ExecutionPolicy::Auto and Cpu. Resident PCG validation policies are rejected rather than silently reused for FGMRES. Adaptive Hybrid/Schwarz/coarse strengthening is not active on the GeneralSquare path yet; the prepared FGMRES workspace is reused across right-hand sides.

To avoid retaining an unused full PCG scratch set, a GeneralSquare prepared context creates a zero-length legacy PCG workspace and allocates the real FgmresWorkspace. SPD and structural execution are unchanged.

0.8-e5 prepared GeneralSquare ILU(0)

GeneralSquare now has an explicit solver-level fixed-preconditioner policy independent from GeneralSquareOptions and its FGMRES restart policy. Jacobi remains the default. Ilu0 is an opt-in prepared preconditioner selected with set_general_square_preconditioner_policy and reused across right-hand sides while the exact matrix structure/value signatures remain unchanged.

Csr32Matrix does not require sorted or unique columns, so ILU(0) owns a private canonical CSR representation. Preparation sorts each row, sums duplicates, removes exact-zero off-diagonal entries, verifies a diagonal entry, and then performs a no-fill numeric LU factorization on that canonical pattern. The caller’s matrix and selected SpMV backend are unchanged. Diagonal factor positions use u32, matching CSR32’s nnz representation.

No row pivoting or fill is introduced. A factor pivot whose magnitude is at or below 1e-12 * row_max is selectively replaced by a same-sign row-relative floor (positive for exact zero). E5 pivot studies found this smallest tested floor sufficient to rescue exact factor-pivot breakdowns, including row scales spanning 1e-8 through 1e8, while applying zero adjustments to the validated healthy hard-B family. The prepared context exposes adjusted-pivot count so callers can detect whether stabilization occurred.

E5 screening also showed that once ILU(0) is active, large Arnoldi restart dimensions can become counterproductive: restart 3 was the fastest and lowest-memory tested setting across the synthetic hard A/B/C families. Nevertheless HyBIT does not couple these controls implicitly. Selecting ILU(0) does not change GeneralSquareOptions::restart; applications must opt into a smaller restart explicitly. Jacobi + fixed restart 30 remains the default GeneralSquare behavior.

0.8-e4b budget-aware GeneralSquare restart control

The opt-in BudgetAware GeneralSquare policy adapts FGMRES restart at exact restart boundaries while remaining inside one Krylov invocation. The Krylov layer exposes restart-boundary progress and accepts the next restart dimension from a controller. This avoids re-entering FGMRES every cycle, preserves the global flexible-preconditioner iteration index, and reuses the one prepared maximum-capacity FgmresWorkspace.

The initial controller is deliberately small and deterministic. It observes four restart cycles at a time, estimates the per-iteration logarithmic residual decay rate, and projects the iterations still required to reach the configured tolerance. Ordinary restart growth requires both sustained projected-budget pressure and a non-improving decay-rate trend. A stronger emergency pressure can override the trend gate. Two consecutive decision windows are required before growth. Restart doubles up to max_restart.

The fixed numeric controller constants are an internal 0.8 checkpoint, not a new tuning surface: ordinary pressure 1.10, emergency pressure 1.40, four cycles per window, and two qualifying windows. E4b studies found that exposing or selecting a pressure threshold alone was not robust across phase-shifted and anisotropic nonsymmetric test families. The trend gate prevented premature growth on cases whose small-restart convergence rate was still improving.

The default remains Fixed; Escalating remains the deterministic E4a alternative. Problems that remain far from tolerance even after restart growth are not treated as restart-controller failures and are expected to feed later nonsymmetric preconditioner escalation work.

0.8-e4a staged GeneralSquare restart escalation

GeneralSquare FGMRES now has an explicit restart-policy layer. Fixed preserves the E2 behavior. The opt-in Escalating policy carries the current solution between FGMRES stages, doubles the restart dimension from restart up to max_restart, and assigns escalation_stage_iterations to each non-final stage. Once the maximum restart is reached, that final stage receives the entire remaining SolverOptions::max_iterations budget.

The prepared context allocates one FgmresWorkspace at the maximum required restart capacity. Lower-restart stages reuse that same Arnoldi/V/Z storage; the workspace does not reallocate at stage boundaries. Matrix analysis, Jacobi preconditioner state, and selected sparse backend are also reused.

This checkpoint is scheduled escalation, not residual-triggered adaptation. E3 experiments showed that small restart dimensions minimize orthogonalization cost on easier nonsymmetric problems, while a harder nonnormal problem benefited from progressively larger restart spaces. The concrete experimental stage lengths motivate this mechanism but are not assumed to be universal.

The default remains GeneralSquareRestartPolicy::Fixed with restart 30, so E4a does not silently change E2 solver behavior.

Future ABTM topology layer

The next ABTM phase is broader than an alternative SpMV storage format. HyBIT will treat bitmap-oriented structure as a symbolic/topology layer that can feed different prepared numerical layouts. Boolean structural operations may prune numerical work for region extraction, masked operators, ILU(0) intersections, coarse construction, and later partition/halo work.

The same logical topology may prepare differently for scalar CPU, block-FEM, Rayon, and GPU execution. ABTM does not claim to eliminate the irregular x[j] gather of ordinary SpMV or guarantee GPU coalescing merely by storing bitmap metadata. Ordering, occupancy, blocking, and target-specific layout are measured parts of preparation.

See ABTM topology algebra and metadata-first execution.

Generic one-shot and prepared path

HybitSolver::solve_csr32 remains available and internally follows analyze -> prepare -> solve. For repeated right-hand sides, the prepared path reuses matrix-dependent state and Krylov workspace.

1. Analyze CSR32 structure, SPD baseline requirements, backend policy, and exact structure/value signatures.
2. Prepare Jacobi and a reusable PCG workspace; ABTM is built eagerly only when required by the selected backend and otherwise lazily on hybrid escalation.
3. Solve #1 with the generic adaptive PCG path.
4. If progress is poor, identify difficult DOFs, form bounded local regions, build local Cholesky corrections, and restart PCG with the strengthened preconditioner.
5. Cache learned local factors for later right-hand sides against the unchanged matrix.

ABTM remains internal; callers provide ordinary CSR32 matrices.

Structural path

solve_structural_csr32 / prepare_structural_csr32 accept 3-D node coordinates matching the CSR DOF ordering and build a six-rigid-body-mode two-level preconditioner. At the r25 checkpoint Structural Auto resolves four independent choices:

• Aggregation: Graph-connected rigid-body aggregates by default, with conservative fallback to contiguous RCM-order aggregates when the graph coarse space cannot be formed safely.
• Fine-grid operator: serial CSR or Rayon-parallel CSR through StructuralSpmvPolicy.
• Preconditioner kernels: serial or parallel 3x3 block-Jacobi, restriction, and prolongation through StructuralPreconditionerPolicy; the packed dense coarse triangular solve remains serial.
• PCG dense vectors: serial or parallel/fused dot/norm/update kernels through StructuralPcgVectorPolicy.

The production structural pipeline is therefore:

CSR32 stiffness + node coordinates
        |
        v
analyze
        |
        v
prepare structural state
  |-- graph/contiguous aggregation
  |-- six rigid-body modes per aggregate
  |-- E = Z^T A Z
  |-- packed Cholesky(E)
  |-- 3x3 block-Jacobi factors
  |-- cached aggregate->node index for parallel restriction
  |-- reusable PCG workspace
        |
        v
solve #1..N
  |-- selected serial/parallel CSR SpMV
  |-- additive rigid-body two-level preconditioner
  |-- selected serial/parallel-fused PCG vector kernels
  +-- reuse prepared factors/index/workspace

Allocation policy

Repeated generic and structural PCG execution reuse preallocated work vectors. Local Schwarz factors own fixed scratch, and the structural parallel restriction index is created during preparation and reused by solve-many. No per-iteration Vec allocation is intended in the production Krylov loops.

Reuse safety

Prepared contexts validate CSR structure and exact floating-point coefficient bits before reuse. This remains intentionally conservative: changing matrix values requires a new prepared context even when the sparsity pattern is unchanged.

0.6 r25 checkpoint

On the development 358065-DOF / 28239653-nnz L-angle physical-load case, 8 Rayon workers resolved Structural Auto to Graph aggregation + Parallel SpMV + Parallel preconditioner + Parallel/fused PCG vectors. The run used a 1398-dimensional coarse space, converged in 220 iterations, and independently verified a relative residual of 9.378557e-9; solve time was 1.726 s and analysis+prepare+solve was 2.797 s. This is a machine-specific regression reference rather than a general performance guarantee.

ABTM topology algebra and metadata-first execution

Development architecture direction for develop/0.8.0. This document describes a planned execution model; it does not claim that ABTM is already the production default or universally faster than CSR.

Purpose

HyBIT began from the ABTM idea: encode sparse topology with bitmap-oriented metadata while keeping numerical values in compact streams. The current solver has since gained prepared execution, hybrid/local correction, algebraic coarse spaces, resident CPU execution, FGMRES, ILU(0), and ordering experiments.

That broader solver architecture changes how ABTM should be positioned.

ABTM is not treated as a universal replacement for CSR. Its primary role is a symbolic/topology representation from which CPU, GPU, block, local, and distributed prepared layouts can be derived.

The central execution principle is:

Metadata-first numerical execution: determine the set of useful numerical operations from topology metadata before loading and combining numerical values.

This separates structural algebra from numerical algebra.

Structural and numerical layers

The logical model is:

Sparse topology              Numerical values
      |                            |
      |-- adjacency                |
      |-- masks                    |
      |-- permutation              |
      |-- regions                  |
      |-- halo                     |
      |                            |
      +-------------+--------------+
                    |
                    v
              prepared layout
                    |
       +------------+-------------+
       |            |             |
     CSR32       CPU ABTM      GPU ABTM

SparseTopology is the reusable structural object. A concrete prepared backend may use CSR32, bitmap chunks, block-ABTM, or another representation without changing the mathematical problem class.

Ordering belongs before layout preparation:

input structure
      |
      v
symbolic graph
      |
      v
ordering / partitioning
      |
      v
ABTM or CSR prepared layout

F1/F2 ordering measurements already show why this matters: the same matrix can have materially different local structure after permutation, which affects both ILU(0) quality and the locality/occupancy of a future bitmap layout.

What ABTM removes, and what it does not

Bitmap topology can reduce or eliminate repeated loading and comparison of explicit column-index metadata. A chunk can represent many possible structural positions with one integer mask, and set bits can be enumerated using standard integer operations.

ABTM does not remove every irregular access. In an SpMV,

y[i] += a_ij * x[j]

the x[j] gather remains unless ordering, tiling, blocking, or local caching also improves vector locality. Therefore the architecture does not claim that ABTM guarantees fully streaming SpMV, branch-free execution, GPU coalescing, or absence of divergence.

The narrower and testable claim is:

• matrix-side topology metadata can be compact and bit-oriented;
• packed numerical values can be streamed in topology order;
• structural filtering can happen before value loads;
• ordering/tiling can improve vector-side locality;
• specialized prepared layouts can map the same logical topology to different CPU and GPU execution strategies.

Topology algebra

The topology layer should expose a small set of operations with clear semantics:

AND       intersection / common support
OR        union / graph expansion
AND-NOT   exclusion / unvisited frontier
XOR       structural difference
popcount  cardinality / work estimate
rank      packed-value offset before a bit
select    bit position from packed ordinal

On Rust scalar backends, use stable integer methods such as u64::trailing_zeros() and u64::count_ones(). Architecture-specific SIMD or std::arch should be introduced only when measurement justifies it.

A statement such as “64 relationships in one clock” is intentionally avoided. A 64-bit logical instruction can process up to 64 Boolean topology states at once, but actual latency and throughput depend on the target processor, dependencies, memory behavior, and compiler code generation.

G1 scalar topology contract

The first implementation checkpoint uses a structural sparse-of-bitmaps object, AbtmTopology, independent of numerical values.

The G1 contract is deliberately structural:

• every stored CSR column position becomes a topology bit, including an explicitly stored numerical zero;
• duplicate stored columns collapse to one topology bit;
• topology construction does not inspect the CSR value array;
• each row stores only non-empty 64-column words;
• row pointers and word indices are u32, while each topology mask is u64;
• the structure-of-arrays payload is therefore 12 bytes per non-empty word plus the row-pointer array;
• AND, OR, AND-NOT, and XOR preserve canonical increasing word order;
• popcount is exact structural cardinality;
• rank counts set bits strictly below a bit/column;
• select maps a zero-based packed ordinal back to a bit/column.

This structural definition is intentionally distinct from the current AbtmMatrix numerical packing, which may discard explicit numerical zeros or duplicate contributions that cancel. Later prepared numerical layouts must state explicitly which topology they consume rather than silently conflating stored structure with numerical activity.

G1 is scalar reference semantics. SIMD, Rayon, GPU preparation, packed numeric updates, and backend promotion are later checkpoints and must not change these Boolean invariants.

Metadata-first product pruning

Consider a sparse dot product

a^T b = sum_k a_k b_k

with chunk masks M_a and M_b. Before touching values:

M_active = M_a AND M_b

If M_active == 0, the complete numerical contribution of that chunk is zero. Otherwise popcount(M_active) gives the exact number of scalar products that remain.

For packed values, the bitmap position must map to its packed ordinal. For a bit position p, conceptually:

rank(mask, p) = popcount(mask AND bits_below(p))

Together with the chunk’’s value offset, rank maps a structural bit directly to its packed numerical value without an index search.

The general execution pattern is:

metadata intersection
        |
        +-- empty -> skip numerical work
        |
        v
     popcount
        |
     work size
        |
        v
   rank / select
        |
        v
 packed value loads
        |
        v
 numerical kernel

This is the key distinction between ABTM as a compressed storage format and ABTM as an execution architecture.

Candidate kernels

Sparse dot and masked reductions

These are the simplest demonstrations of metadata-first pruning. The bitmap intersection defines exactly which products are required.

Masked and restricted SpMV

Ordinary SpMV has only the matrix support mask and therefore cannot use a second mask to remove arbitrary products. Restricted operators do have another mask:

effective = row_topology AND active_dof_mask

This applies naturally to local residuals, subdomain operators, boundary-only operations, constrained/free DOF masks, and selected hard regions.

Sparse matrix products

For

C_ij = sum_k A_ik B_kj

a row-support mask of A and column-support mask of B can determine whether the dot product is structurally empty before any numerical values are loaded. A prepared matrix-product backend may therefore retain a transpose/column topology in addition to row topology. Numerical values need not necessarily be duplicated merely because topology is available in both orientations.

Region growth

Hard-region expansion can use:

next = neighbors(frontier) AND-NOT visited

where local adjacency masks are OR-combined before the exclusion step. This is a topology operation and should not require numerical matrix values.

Local submatrix extraction

For a region mask R and matrix row topology A_i:

local_i = A_i AND R

popcount(local_i) gives the local row nnz before allocation; rank/select then extract only the required packed values. This directly supports bounded local direct factors used by HyBIT’’s escalation path.

Region overlap and multiplicity

Pairwise overlap is a bitmap intersection, but Schwarz weight multiplicity is not merely a pairwise Boolean question when more than two regions share a DOF. The implementation should build explicit multiplicity state, for example a small integer counter per participating DOF or a bit-sliced counter, from the region masks. Numerical matrix values are still unnecessary.

ILU(0) symbolic/numeric intersection

For canonical no-fill ILU(0), an update can be restricted to structural intersections between the current row and the relevant upper-factor row. A future ABTM ILU(0) experiment should compare:

CSR search / merge

against

bitmap intersection -> rank -> packed update

rather than assuming the bitmap route is faster. This is especially relevant to the current GeneralSquare work because ILU(0) is already a validated prepared preconditioner.

Coarse/Galerkin construction

The support of coarse couplings in P^T A P can be screened through topology before accumulating numerical products. ABTM may therefore accelerate symbolic coarse construction even when the final coarse matrix is stored in a different format.

MPI halo extraction

For a partition, topology can distinguish local/internal and remote/boundary dependencies before numerical execution. Bitwise union/intersection is a candidate for deriving halo DOFs and communication schedules; overlap of nonblocking communication with internal computation remains a later distributed-memory concern.

Adaptive and block ABTM

A fixed bitmap representation is not expected to dominate every sparsity pattern. Chunk occupancy determines whether bitmap metadata is cheaper than explicit indices.

A prepared implementation should be able to choose among representations such as:

enum TileKind {
    Dense,
    Bitmap,
    Sparse,
    Block3,
    Block6,
}

This is an architectural direction, not a frozen public API.

For FEM, node-level topology is especially important. A single node-adjacency bit can guard an entire dense 3 x 3 or 6 x 6 coupling block. The metadata cost is then amortized over 9 or 36 scalar coefficients and the numerical kernel can use fixed-size block loads/FMA sequences.

CPU and GPU preparation

A logical ABTM topology should not force one physical descriptor layout on all targets.

CPU preparation may favor 64-bit chunks, compact descriptors, rank/popcount, and Rayon row/tile partitioning.

GPU preparation may instead favor warp-sized work units, bounded bitmap-word counts, cooperative decode, local/shared caching of vector tiles, and device-specific descriptor packing. A future CubeCL backend should therefore consume a GPU-prepared ABTM representation rather than copying the CPU layout verbatim.

Coalesced access and low divergence are optimization targets to measure, not properties guaranteed by the logical bitmap representation.

Dynamic topology

When nonlinear analysis or contact changes sparsity, structural edits should be buffered separately from packed numerical storage. A rebuild can merge topology changes first and repack values in one preparation step instead of repeatedly reallocating the complete numerical stream.

The prepared-state API should distinguish:

topology unchanged, values changed
topology changed

so later symbolic/topological reuse can be separated from numeric refactorization.

Measurements required

ABTM evaluation should not be reduced to “SpMV faster or slower than CSR”. Each experiment should record at least:

• chunk/tile occupancy;
• metadata bytes per nnz;
• packed-value bytes;
• bitmap operations performed;
• candidate numerical products;
• executed numerical products;
• products skipped by metadata;
• rank/select cost;
• setup/preparation cost;
• solve/kernel wall time;
• cache/locality counters where practical.

Define a pruning ratio:

pruning_ratio =
    1 - executed_numerical_products / candidate_numerical_products

The useful break-even condition is qualitative but fundamental:

metadata cost < numerical and memory work avoided

A metadata intersection that avoids one cheap FMA may lose. The same intersection that discards many scalar products, a 3 x 3/6 x 6 block, a local-factor update, or a remote communication dependency can be highly profitable.

Development sequence

After the current GeneralSquare robustness checkpoints, the ABTM track should proceed experimentally:

1. G1 — topology algebra: bitmap chunks, AND/OR/AND-NOT/XOR, popcount, rank/select, invariants, and scalar reference tests.
2. G2 — metadata-first kernels: sparse dot/support intersection and quantitative pruning metrics.
3. G3 — HyBIT region operations: region growth, overlap/multiplicity, and local submatrix extraction versus the current graph/index path.
4. G4 — ABTM ILU(0): symbolic intersection and packed numeric updates versus canonical CSR ILU(0).
5. G5 — block ABTM (validated): 3x3/6x6 FEM node topology and explicit fixed-size dense block-CSR numerical execution.
6. G6 — prepared restricted SpMV (validated): ordinary scalar ABTM remains diagnostic; fixed AM and graph-local RA*R restrictions materialize compact CSR, with preparation policy kept explicit and workload-specific.
7. G7 — Rayon prepared execution (validated): explicit full-Rayon and task-limited prepared CSR execution, serial gather, and no automatic hardware-specific size selector.
8. G8 — distributed topology (next): partition/halo extraction and MPI-host integration before distributed Krylov scheduling.
9. G9 — GPU prepared execution: CubeCL-oriented descriptors and resident execution after the CubeCL API is sufficiently stable.

No automatic backend promotion should occur until the relevant workload class shows repeatable benefit including preparation cost.

Architectural invariant

The long-term invariant is:

Boolean topology algebra determines where work can exist; numerical algebra is invoked only for the surviving work, using a representation prepared for the current execution target.

This makes ABTM a symbolic computation layer for HyBIT rather than merely another sparse matrix file format.

G1 validation evidence and decision

G1 establishes scalar reference semantics for structural bitmap topology and validates the representation on the same ten-matrix corpus used by the recent GeneralSquare studies.

The measured AbtmTopology physical form is a sparse-of-64-bitmaps structure-of-arrays:

row pointer  : u32 per row boundary
word index   : u32 per non-empty word
word mask    : u64 per non-empty word

so one non-empty bitmap word costs 12 bytes before row-pointer amortization.

Ten-matrix metadata/occupancy evidence

Eight of ten matrices use less topology metadata than CSR32 row-pointer plus column-index metadata. sherman5 is approximately break-even and thermal1 is materially worse because most 64-column words contain very few structural entries.

The result confirms that the logical bitmap algebra is broadly useful, but a single physical bitmap encoding is not an appropriate universal prepared layout.

Occupancy and adaptive representation

The observed occupancy distribution is highly matrix dependent.

• thermal1 has median occupancy 1 and about 83% of words at three entries or fewer.
• sherman5 has median occupancy 3 and about 81% of words at three entries or fewer.
• raefsky3 has median occupancy 24.
• nd3k has median occupancy 17 and about 11% of words at 40 entries or more.
• cant, s3dkq4m2, and x104 occupy an intermediate regime.

For a conceptual sparse word representation

word_index: u32
offsets:    k x u8

the metadata model is roughly 4 + k bytes, compared with 12 bytes for the bitmap word. This gives a natural conceptual crossover near k = 8, consistent with the existing ABTM sparse/bitmap design direction. Concrete prepared layouts must still account for offsets, alignment, descriptor storage, and execution cost rather than freezing this arithmetic as a universal threshold.

Rank/select interpretation

G1a deliberately checked both chunk-local and row-wide convenience rank/select. On dense rows this made the row-wide scan visible, most strongly on nd3k.

G1b isolates the operation required by packed numerical streams: chunk-local word rank/select. Across the ten real matrices it measured approximately 3.5–9.3 ns per structural nonzero on the validation machine. Therefore future metadata-first kernels should carry the matching word descriptor and use word-local rank/select; row-wide rank/select remains a correctness/convenience API rather than the inner numerical addressing path.

Boolean merge validation

G1b also prepared two deterministic approximately 75%-dense structural subsets of each source topology and exercised AND, OR, AND-NOT, and XOR through the general merge path.

The measured intersection ratios stay close to the independent expectation

0.75 * 0.75 = 0.5625

and union ratios stay close to

1 - 0.25 * 0.25 = 0.9375

across all ten matrices. This validates missing-word and partial-mask merge semantics in addition to the unit-test identities and G1a self-operations.

G1 decision

G1 is accepted as the scalar logical topology layer.

The retained invariants are:

• numerical values are not part of topology;
• explicitly stored structural positions remain topology even when their numerical value is zero;
• duplicate structural columns collapse to one bit;
• row words remain canonical and strictly increasing;
• AND / OR / AND-NOT / XOR operate only on topology;
• popcount is exact structural cardinality;
• word-local rank/select is the packed-value addressing primitive.

G1 does not promote one bitmap physical layout as universal. Later prepared execution may choose Sparse, Bitmap, Dense, block, CPU-specific, or GPU-specific representations from the same logical topology.

G2-G7 are now validated. The next checkpoint is G8: distributed topology / MPI-host integration.

G2a scalar metadata-first sparse-dot experiment

G2 begins with a deliberately narrow numerical kernel: one sparse matrix row dotted with a vector whose structural support is represented by a DofMask.

The scalar prepared reference stores one numerical value per G1 structural position in the same packed topology order. For each 64-column topology word:

active = matrix_word_mask AND vector_support_word

If active == 0, no numerical value from that matrix word is loaded. Otherwise word-local rank maps each active bit directly to its packed matrix value.

For this experiment the work metrics are defined precisely:

candidate_products = structural matrix entries before support filtering
executed_products  = popcount(matrix topology AND vector support)
skipped_products   = candidate_products - executed_products
pruning_ratio      = skipped_products / candidate_products

This makes G2a a sparse-dot/support-intersection experiment rather than a claim about a production ABTM SpMV implementation. Full masked/restricted SpMV policy, parallel execution, and backend selection remain later checkpoints.

The 100%-active case measures metadata/rank overhead when no numerical work can be pruned. Lower support densities measure whether avoided value loads and multiplications can amortize that overhead.

G2c dual row/column topology

Matrix-versus-matrix sparse products require metadata on both operands. For

C_ij = sum_k A_ik B_kj

the structural intersection is between a row support of A and a column support of B.

G2c introduces AbtmDualTopology, which retains:

row topology       original matrix structure
column topology    structural transpose topology

while leaving numerical values unduplicated. A future numerical prepared backend can therefore decide independently whether values also need a transpose-oriented stream.

The G2c benchmark first isolates structural work. It compares a conventional sorted explicit-index intersection with bitmap-word intersection for deterministic row/column pairs. The key metrics are metadata footprint, metadata comparisons, mask AND count, empty-dot ratio, exact overlap products, and scalar intersection wall time.

This avoids conflating topology benefit with a particular numerical-value layout before the dual-metadata economics are known.

G2d numerical values for dual topology

G2c establishes that row/column bitmap metadata can greatly reduce structural intersection work on high-occupancy matrices, while low-occupancy cases remain neutral or unfavorable.

G2d asks the next architectural question: must numerical values also be duplicated for the column orientation?

Two bitmap variants are compared:

mapped-column values
    one packed row-value stream
    + u32 source-value index per column-topology entry

duplicated-column values
    packed row-value stream
    + packed column-value stream

Both use the same dual row/column topology. A conventional explicit row/column-index merge with duplicated values is retained as the numerical baseline.

This separates three costs that must not be conflated:

topology compression
column-value addressing indirection
column-value duplication

The result will determine whether later G2 matrix-product kernels should keep a single numerical copy, duplicate values for hot transpose-oriented work, or select between the two as a prepared policy.

G2e adaptive numerical dual layout

G2d shows that one fixed bitmap numerical addressing strategy is insufficient: high-occupancy matrices can win, but low- and medium-occupancy matrices pay too much rank/value-addressing overhead.

G2e therefore returns to the physical-layout rule established by G1. Both row and column orientations use the existing adaptive tile classification:

Sparse <= 8 entries / 64-column word
Bitmap intermediate occupancy
Dense  >= 40 entries / 64-column word

Execution is also adaptive. Sparse/Sparse tiles use a local ordered set-bit merge and avoid rank. Dense tiles use direct bit offsets. Bitmap tiles retain mask intersection and packed rank addressing.

This checkpoint tests whether representation and execution policy must both be occupancy-aware. The thresholds are intentionally left unchanged until the mechanism is measured on the corpus.

G2f physically packed adaptive metadata

G2e demonstrates that occupancy-aware execution materially improves every tested matrix compared with fixed bitmap numerical addressing. However, the existing TileDesc is still a fixed 16-byte descriptor. In particular, a Sparse tile still carries a full 64-bit mask and value offset.

G2f separates the representation question from the execution question. Sparse tiles use variable-length byte offsets, while Bitmap and Dense tiles retain 64-bit masks. Row-local payload and value pointers make each row directly addressable without storing a value offset in every tile.

The benchmark compares the resulting dual numerical layout directly with explicit CSR plus transpose-CSR storage and numerical row/column dot time. Threshold tuning is intentionally deferred until the packed representation itself is validated.

G2g typed compact descriptors

G2f establishes that compact Sparse storage is possible but also shows that a byte-oriented variable payload should not be decoded directly in the numerical hot loop.

G2g therefore evaluates an intermediate design:

8-byte typed descriptor + typed auxiliary streams

Sparse offsets remain one byte each and Bitmap/Dense masks remain 64-bit, but the descriptor directly identifies the relevant typed side stream. Value offsets are advanced sequentially within each row, preserving compactness without per-tile value offsets.

This checkpoint is intended to find the practical middle ground between G2e’s fast fixed 16-byte descriptors and G2f’s compact but expensive variable decoder.

G2 closeout

G2 is closed for the 0.8 development checkpoint. The full evidence and ten-matrix G2e corpus are recorded in ABTM_G2_CLOSEOUT.md.

The selected direction is dual structural topology plus adaptive Sparse/Bitmap/Dense CPU execution, while retaining explicit CSR/CSC-like fallback. G2 does not promote a universal bitmap layout or an automatic ABTM/explicit threshold.

G3 is the next checkpoint and moves to region growth, overlap/multiplicity, and local submatrix extraction.

G3a region growth

G3 begins with structural region growth. For a seed set R_0, G3a defines undirected adjacency from the matrix structure as

G = support(A) union support(A^T)

and grows only the previous hop’s frontier:

F_0 = R_0
N_k = neighbors(F_k)
F_{k+1} = N_k \ R_k
R_{k+1} = R_k union F_{k+1}

Using AbtmDualTopology avoids constructing numerical transpose values for this operation. Row topology supplies outgoing neighbors and column topology supplies incoming neighbors. The result is deterministic and purely structural.

G3a intentionally does not yet assign multiplicity or extract local matrices; those follow after region-growth semantics and cost are validated.

G3b overlap and multiplicity

For a family of grown regions R_i, G3b defines the node multiplicity

m(v) = sum_i [v in R_i].

This gives three directly useful domain-decomposition quantities:

• covered nodes: m(v) >= 1;
• overlap nodes: m(v) >= 2;
• multiplicity-weighted overlap work.

The implementation retains the full per-node multiplicity map and verifies the identity

sum_{i<j} |R_i intersect R_j| = sum_v m(v)(m(v)-1)/2.

This checkpoint establishes overlap semantics before G3 local submatrix extraction. It does not yet define partition ownership or weighting policy.

G3c structural local-submatrix extraction

For a region R, G3c constructs the CSR-like structural pattern of A[R,R]. Rows and columns use deterministic ascending-global-index local numbering.

ABTM performs this as a metadata-first operation:

for global row r in R:
    for topology word W in row(r):
        kept = W.mask AND R.word(W.index)
        enumerate only kept bits

The extracted pattern stores local row pointers and local column indices plus the ordered global-node list. A temporary dense global-to-local map is currently used and its scratch bytes are reported explicitly.

This checkpoint deliberately stops at structural extraction. Numerical value gathering and local numeric factor preparation remain separate so topology benefits are not hidden by a value-addressing policy chosen too early.

G3d prepared local numeric refresh

G3d separates stable symbolic structure from changing numerical coefficients.

For each region, preparation stores:

local pattern:
    global nodes
    row_ptr
    col_idx

numeric address plan:
    local structural entry -> CSR source value position(s)

If the input CSR is canonical, each local structural entry has one source position and refresh becomes a direct indexed gather. If duplicate CSR entries exist, the plan stores all source positions and sums them without changing the local symbolic pattern.

This design targets repeated refactorization/nonlinear/time-stepping use cases: region growth and local structure are paid once while numerical values can be refreshed many times.

G3 closeout decision

G3 confirms the topology/numerical boundary introduced in G2.

ABTM dual topology is retained for region discovery, frontier set algebra, overlap metadata, and symbolic pruning. Local CSR-like structures are materialized when downstream kernels need explicit row/column indexing. Prepared source-address maps are appropriate when numerical coefficients change while the local symbolic structure remains fixed.

The G3d refresh result is a symbolic-reuse result, not a claim that ABTM bitmap value storage is universally faster than prepared CSR. See ABTM_G3_CLOSEOUT.md for the corpus and interpretation limits.

G4a ILU(0) symbolic intersection

G4 begins by isolating ILU(0) symbolic update lookup from numerical factorization.

For each lower structural entry (i,j), canonical CSR ILU(0) probes the upper structure of row j and binary-searches candidates in the remainder of row i. The G4a experiment replaces those per-entry probes with word-level intersection of the same canonicalized pattern.

The benchmark does not use raw topology directly when explicit zeros or duplicate cancellation would differ from canonical ILU(0). It first applies the same canonical structural rule as the production preconditioner, then builds the ABTM topology.

No production ILU(0) routing changes are made in G4a.

G4b ILU(0) numeric intersection

G4b adds a compact row-word prefix to turn topology intersection bits into canonical CSR numerical positions. This avoids a per-successful-update address plan, whose storage could grow with the number of ILU(0) products rather than with matrix topology.

The experiment keeps numerical factors in conventional contiguous f64 CSR order. ABTM supplies symbolic intersection and rank-within-word addressing; it does not replace the numerical factor layout.

G4c ILU(0) rank-LUT addressing

G4c tests a direct rank lookup for numerical addressing. Each nonempty topology word receives a 64-byte table mapping a set bit position to its ordinal among the word’s structural bits.

This removes two hardware popcount rank operations per executed ILU(0) update, at the cost of 64 bytes per nonempty word. G4c is diagnostic: if the speed gain is small or inconsistent, the extra storage is rejected. If the gain is significant only for dense words, a later adaptive threshold may retain tables only where occupancy justifies them.

G4d adaptive rank-LUT sweep

G4d evaluates selective rank LUT storage by topology-word occupancy. Dense words may receive a direct 64-byte rank table while sparse words retain rank-by-popcount.

The experiment explicitly measures both static selected-word fraction and dynamic LUT-rank coverage. The distinction matters because words are not accessed uniformly during ILU(0): a structurally modest word may still be hot across many lower-pivot intersections.

No production threshold is implied by average word occupancy alone.

G4f production lifetime boundary

The G4c rank LUT is required only while constructing ILU(0). G4f makes that lifetime explicit: topology words, structural prefixes, and u8[64] rank tables are preparation scratch and are released once the canonical L/U values have been produced.

As a result, ABTM factorization can accelerate setup without increasing the persistent ILU(0) factor footprint or changing triangular application.

G5 block-topology validation and production decision

G5 tests the FEM-oriented hypothesis that scalar topology can be coarsened into node-sized 3x3 or 6x6 couplings and that sufficiently full numerical blocks can amortize index/gather overhead.

G5a topology evidence

Across the 11-matrix development corpus, block topology compressed structural metadata substantially:

The B6 topology is especially compact, but low block fill can make dense value storage prohibitively expensive. Therefore block topology and numerical physical layout must remain separate decisions.

G5b-G5d numerical kernel evidence

A generic runtime-sized block kernel lost to scalar CSR on all 11 development matrices. G5 therefore does not promote a generic block traversal.

Fixed-size B3/B6 kernels changed the result. Tail-specialized execution keeps full interior blocks on the fixed kernel and isolates only the final partial row/column block. Representative kernel speedups versus scalar CSR were:

| d3k | B6 | 0.741 | 2.57x | 0.912 | 62.8 | | x104 | B6 | 0.857 | 2.31x | 0.786 | 66.8 | | aefsky3 | B6 | 0.748 | 2.08x | 0.900 | 73.8 | | cant | B3 | 0.930 | 1.88x | 0.755 | 77.6 | | oneS01 | B3 | 0.822 | 1.47x | 0.853 | 116 | | s3dkq4m2 | B6 | 0.549 | 1.35x | 1.224 | 171 |

Low-fill matrices such as hermal1, pache2, cfd1, and enkat25 remain better served by scalar CSR.

G5e-G5f selector evidence

The development selector was deliberately simple:

ext if B3 fill < 0.65: CSR else if B6_fill / B3_fill >= 0.75: B6 else: B3

On the 11-matrix calibration corpus it produced no false-positive block routes. Its one kernel-winner miss was sherman5, where B3 was only about 1% faster and required roughly 2088 SpMV to repay preparation, so CSR remained the sensible finite-horizon choice.

The thresholds were then frozen. Held-out validation on Goodwin_010, G3_circuit, parabolic_fem, hermal2, and inline_1 classified all five correctly with:

ext false_positive_block = 0 false_negative_csr = 0 block_size_miss = 0

inline_1 provides the positive held-out B3 case: B3 fill about 0.99998, kernel speed about 1.78x CSR, storage about 0.70x CSR, and preparation break-even about 73.5 SpMV.

No positive held-out B6 case is available yet. The structural rule is therefore retained as development evidence, not promoted to an automatic production selector.

Explicit G5 production API

G5 exposes the validated numerical representation explicitly:

` ust use hybit::{DenseBlockCsrOperator, DenseBlockSize};

let block = DenseBlockCsrOperator::from_csr32(&matrix, DenseBlockSize::B3)?; `

The prepared numerical layout is conventional compact block CSR:

ext block row_ptr block col_idx fixed dense BxB values

The operator supports B3/B6, partial final row/column blocks, duplicate scalar contribution accumulation, storage/fill diagnostics, and the common LinearOperator interface.

This preserves the ABTM architectural boundary: topology and measurements identify profitable structure, while the hot numerical representation can be a specialized conventional layout. Existing CSR paths and solver defaults are unchanged.

G6 is the next checkpoint: ordinary and masked/restricted SpMV.

G6 validation evidence and decision

G6 evaluated ordinary SpMV, masked column restriction, and graph-local restriction as separate workload classes instead of forcing one ABTM numerical layout onto all three.

G6a ordinary and on-the-fly masked SpMV

The 11-matrix development corpus rejected scalar ABTM as a universal ordinary SpMV replacement. Ordinary ABTM produced no wins and a geometric-mean ABTM-versus-CSR speedup of about 0.43x.

Metadata-first masking still demonstrated useful pruning semantics. Intersecting tile topology with the active-DOF mask before loading numerical values strongly beat a CSR implementation that branches on every stored entry. However, the stronger baseline–a pre-zeroed input vector followed by ordinary CSR SpMV–was harder to beat and did not justify promoting the on-the-fly ABTM masked kernel.

The retained conclusion is that metadata-first pruning is valuable during preparation, but repeated scalar arithmetic should use a layout specialized for the prepared workload.

G6b fixed column restriction A*M

For a fixed active-column mask M, G6b materialized a compact CSR operator that retains global dimensions but stores only active-column entries. Repeated apply therefore consumes the changing global input vector directly without a per-entry mask branch or a per-apply masked-vector refresh.

Development-corpus geometric-mean speedup of prepared execution versus the dynamic full-CSR baseline was approximately:

When an AbtmMatrix already existed, ABTM metadata intersection sometimes accelerated preparation. Development-corpus preparation speedup versus direct CSR scan was geometrically about 1.17x--1.35x through 50% density. The advantage was not universal, so it remains explicit rather than automatic.

G6c graph-local R*A*R

G6c built deterministic graph-local regions and compared repeated scanning of global CSR region rows against an explicitly materialized compact local CSR.

Development-corpus dynamic prepared speedup was:

Unlike the column-mask case, the existing ABTM local numeric-plan builder lost to direct CSR region-row extraction on every measured development case. The local plan performs general topology extraction and source binding work that is unnecessary when a one-shot compact R*A*R operator can be built directly.

G6d held-out validation

The held-out corpus was Goodwin_010, G3_circuit, parabolic_fem, thermal2, and inline_1.

For fixed A*M, prepared dynamic execution won 5/5 held-out matrices at 5%, 10%, 25%, and 50% density and 4/5 at 75%. Geometric-mean speedups were about 4.249x, 2.798x, 2.048x, 1.745x, and 1.229x respectively.

The held-out ABTM-preparation advantage was mixed: geometric-mean ABTM-versus-direct-CSR preparation speedup was about 1.03x, 1.00x, 0.97x, 1.07x, and 0.87x over the same densities. This rules out an unconditional ABTM preparation selector.

For graph-local R*A*R, prepared dynamic execution won 5/5 held-out matrices at every tested region fraction. Geometric-mean speedups at 1%, 5%, 10%, and 25% were about 3.727x, 2.682x, 2.516x, and 2.381x.

Direct CSR local preparation beat the existing ABTM local-plan path 5/5 at every held-out fraction, by geometric means of about 6.84x, 7.82x, 8.09x, and 8.57x.

G6 production boundary

G6 therefore retains this workload-specific boundary:

ordinary SpMV
    -> CSR remains the baseline numerical path

fixed column restriction A*M
    -> compact prepared CSR
    -> direct CSR preparation by default
    -> explicit ABTM preparation only when an ABTM layout already exists

graph-local R*A*R
    -> direct CSR region-row extraction
    -> compact local CSR

repeated arithmetic
    -> conventional compact CSR

The public production surface is explicit:

• PreparedColumnRestrictedCsrOperator
  • from_csr32
  • from_abtm as an explicit opt-in preparation path
• PreparedLocalCsrOperator
  • from_csr32
  • compact ascending local numbering
  • gather_input / gather_input_vec for changing global vectors

No automatic restriction policy, CSR/ABTM preparation selector, or solver backend promotion is introduced by G6.

This reinforces the ABTM architectural role: topology metadata is used where it can eliminate structural work profitably, while repeated numerical execution is free to use a conventional representation better matched to the target workload.

G7 validation evidence and decision

G7 evaluates CPU parallel execution only after G6 has already materialized the repeated numerical workload as compact prepared CSR. The question is therefore not whether ABTM itself should become the parallel numerical storage format, but how prepared A*M and R*A*R operators should expose Rayon execution without hiding hardware-specific policy.

G7a prepared execution crossover

With the default 16-worker Rayon pool, fixed-column A*M was strongly parallel-friendly once the prepared matrix was beyond the smallest test case. Across the 11-matrix development corpus the geometric-mean parallel speedups versus serial prepared CSR were approximately:

The five losses were the very small sherman5 prepared operators. The development envelope was max_loss_nnz=15,330 and min_win_nnz=28,316, which suggests a useful size signal but does not by itself define a portable production threshold.

For graph-local R*A*R, full-pool parallelism was more size-sensitive. At 16 workers the kernel geometric-mean speedups at 1%, 5%, 10%, and 25% region fractions were about 0.298x, 1.178x, 1.814x, and 3.118x. The corresponding changing-global-vector dynamic path, including serial gather_input, measured about 0.331x, 1.112x, 1.580x, and 2.015x.

G7b worker-count sensitivity

A 4/8/16-worker sweep showed that A*M continues to favor wider parallelism as the prepared matrix grows. Across 55 development cases, the fastest measured worker count was 4, 8, and 16 in 8, 22, and 25 cases respectively.

R*A*R behaves differently. Across the 44 development cases, the fastest kernel worker count was 4/8/16 in 20/17/7 cases, and the fastest dynamic gather-plus-SpMV path was 4/8/16 in 24/16/4 cases.

This demonstrates that one global worker count is not a sufficient portable prepared-local execution policy.

G7c task granularity and gather policy

G7c retained the default 16-worker global Rayon pool and varied only the number of contiguous row tasks submitted per operation: 2, 4, 8, or 16.

The result rejected the idea that a small fixed task count can simply reproduce the favorable 4-worker behavior. Sixteen tasks were the fastest full dynamic path in 39 of 44 development cases; 8 tasks won 4 cases, 4 tasks won 1, and 2 tasks won none.

The experiment also isolates global-to-local gathering. Parallel gather was usually slower than the existing serial PreparedLocalCsrOperator::gather_input. Even at 16 tasks, the gather geometric-mean speedup versus serial gather was only about 0.030x, 0.126x, 0.246x, and 0.545x at the tested 1%, 5%, 10%, and 25% region fractions. Therefore G7 retains serial gather.

With serial gather retained, 16-task local SpMV became profitable around the low-tens-of-thousands nnz range on the development corpus, but this observation is treated as evidence rather than API semantics.

G7d held-out selector validation

G7d froze a candidate development selector and tested it on Goodwin_010, G3_circuit, parabolic_fem, thermal2, and inline_1.

For fixed-column A*M, the candidate was:

prepared nnz < 32k  -> serial CSR
otherwise           -> existing full Rayon CSR

This candidate produced no held-out regressions at any tested density. The oracle-capture geometric mean was 1.0 at 5%, 10%, 25%, 50%, and 75% active columns, meaning the selector chose the fastest measured serial/full-Rayon path for every held-out A*M point.

For graph-local R*A*R, the candidate was:

local nnz < 32k        -> serial gather + serial CSR
32k <= nnz < 512k      -> serial gather + capped contiguous-row tasks
local nnz >= 512k      -> serial gather + existing full Rayon CSR

The aggregate result was positive but not strong enough to freeze the policy. At 5%, 10%, and 25% region fractions the candidate had no meaningful held-out regression, but at 1% only 4/5 cases were non-regressions.

Two cases explain why an automatic selector remains deferred:

• G3_circuit, 1% region, 74,239 local nnz: the candidate chunked path ran at about 0.848x the serial baseline.
• inline_1, 1% region, 348,237 local nnz: the candidate chunked path was still about 1.97x faster than serial, but the existing full-Rayon path was much faster, reducing oracle capture to about 0.608.

The held-out evidence therefore supports explicit Rayon capability but not a portable hardware-independent R*A*R nnz selector.

G7 production boundary

G7 exposes execution mechanisms without embedding the benchmark machine’s routing thresholds:

Csr32Matrix
    apply_parallel(...)
    apply_parallel_with_tasks(...)

ParallelCsr32Operator
    apply(...)
    apply_with_tasks(...)

PreparedColumnRestrictedCsrOperator
    apply(...)                       # serial LinearOperator path
    apply_parallel(...)
    apply_parallel_with_tasks(...)

PreparedLocalCsrOperator
    gather_input(...)                # serial
    apply(...)                       # serial LinearOperator path
    apply_parallel(...)
    apply_parallel_with_tasks(...)

target_tasks bounds contiguous row-task granularity while retaining the process-wide Rayon pool. It does not create a private pool and does not promise that exactly that many workers run concurrently.

No automatic matrix-size threshold, hardware-specific worker policy, or prepared-operator route is promoted in G7. The caller remains responsible for selecting serial, full-Rayon, or task-limited execution from application and hardware context.

This closes the CPU prepared-execution sequence while preserving the broader ABTM architectural rule: topology discovers and prepares the workload; the numerical representation and execution policy remain target-specific and evidence-driven.

ABTM G2 closeout: metadata-first product pruning

Status: validated / closed for the 0.8 development checkpoint

G2 asked whether ABTM topology can eliminate numerical work before value access, whether row/column topology should be prepared in both orientations, and which numerical physical layout is appropriate for CPU sparse row-by-column products.

The answer is not one universal storage format. The evidence supports a metadata-first topology layer, adaptive numerical execution, and an explicit fallback for structures where ABTM does not pay for itself.

Experiment sequence

G2e ten-matrix numerical corpus

All timings are median scalar CPU row-by-column sparse-dot measurements with eight deterministic pairs per row and five repeats. time ratio and storage ratio are relative to an explicit CSR + transpose-CSR baseline.

The corpus does not justify selecting only from 64-column word occupancy. boneS01 is the counterexample: about 98.6% of its tiles are Sparse, yet the adaptive ABTM path is about 21.8% faster than explicit and uses about 10.7% less storage. Row structural length and the amount of row/column merge work are therefore material predictors in addition to local word occupancy.

Within this corpus, average row length separates the clear outcomes: matrices at or below about 28 nnz/row are explicit or near parity, while all measured matrices at or above about 43 nnz/row favor G2e. The unmeasured region between those groups is too large to promote a production threshold. An avg_nnz_per_row ~= 32 rule is useful only as an experimental screening heuristic, not a stable API or automatic policy.

Numerical agreement

The adaptive and explicit kernels accumulate the same sparse products in different orders. The initial six-matrix corpus had scaled errors from exact agreement through approximately 1.5e-13. The larger held-out matrices exposed more cancellation:

The performance-harness validation tolerance is therefore 1e-8. This is a benchmark-equivalence gate, not a solver convergence criterion. Any production promotion must additionally validate operator-level relative error and end-to-end residual behavior.

Architecture decision

G2 closes with the following CPU direction:

logical/symbolic layer:
    dual row/column topology
    metadata-first intersection and pruning

hot numerical layer:
    fixed direct descriptor
    adaptive Sparse / Bitmap / Dense execution
    prepared row values
    duplicated prepared column values when transpose-oriented reuse justifies it

fallback:
    explicit CSR / CSC-like path

The following are explicitly not promoted:

• universal bitmap-only numerical storage;
• source-index mapped single-copy column values as the hot default;
• byte-decoded variable payloads in the numerical hot loop;
• typed compact 8-byte descriptors as a replacement for the faster G2e layout;
• a production automatic ABTM/explicit threshold derived from this corpus.

Preparation cost

G2e prepares both row and column orientations. On the larger corpus matrices this is on the order of tens to a few hundreds of milliseconds per orientation. That cost must be amortized by repeated row/column work. A production selector therefore needs a solve/work horizon or an application-provided reuse signal; kernel time alone is not sufficient.

Next checkpoint

G3 moves from generic sparse-dot experiments to ABTM’s intended topology-heavy operations:

1. region growth;
2. overlap and multiplicity;
3. local submatrix extraction;
4. preparation/reuse metrics for those operations.

G3 should keep CSR as the numerical reference/fallback and should not assume that the G2e physical layout must be reused unchanged for every region kernel.

HyBIT 0.6 real FEM benchmark design

Goal

Measure solver behavior on assembled, constrained, real SPD FEM systems using a portable input format and a reproducible baseline. This benchmark is intended to answer whether adaptive local direct correction improves convergence and wall-clock behavior on application matrices, not only on synthetic regressions.

Input

fem_bench reads Matrix Market coordinate files:

• field: real or integer;
• symmetry: general or symmetric;
• indices: standard Matrix Market 1-based indices;
• output storage: zero-based CSR32 with sorted columns and duplicate entries summed.

The current automatic solver requires an SPD system. Essential boundary conditions therefore need to be applied before the matrix is exported.

Fair comparison

The benchmark uses the same matrix, zero initial guess, tolerance and iteration limit for both paths:

1. plain CSR32 Jacobi-PCG via hybit::pcg;
2. HyBIT Auto via HybitSolver::solve_csr32.

The plain path is not routed through HyBIT’s probe/restart policy; it is one ordinary PCG recurrence with a fixed Jacobi preconditioner.

RHS modes

By default the benchmark constructs x_exact = 1 and b = A*x_exact. This provides a known solution and isolates the linear solver from load-vector I/O. A physical RHS may instead be supplied with --rhs FILE; the file is a whitespace-separated list of f64 values.

Reported quantities

• dimensions, nnz, average/max nnz per row;
• CSR storage and metadata bytes;
• load time;
• plain Jacobi setup and solve time;
• HyBIT analysis, prepare, probe, diagnostics, local factor and solve time;
• iteration counts;
• solver-reported and independently recomputed relative residuals;
• hard DOFs, number and maximum size of local regions;
• local-factor and Krylov-workspace memory;
• known-solution relative error when using generated RHS.

Interpretation

Do not infer a general speedup from a single run. Large problems should be run multiple times in separate processes, with CPU frequency/power policy and background load controlled. The first 0.6 experiments should focus on correctness, convergence behavior and whether selected regions are numerically meaningful before tuning for throughput.

Structural block-Jacobi diagnostic

For vector-valued structural FEM matrices whose free displacement DOFs remain contiguous per node, HyBIT includes a separate block-Jacobi benchmark. This intentionally stays outside HybitSolver Auto policy until it has been measured on real matrices.

For the L-angle export, use three displacement DOFs per node:

.\benchmarks\scripts\bench-fem-block.ps1 D:\path\to\L-angle-K.mtx -BlockSize 3

The benchmark uses the same known-solution construction b=A*1, factors each contiguous diagonal block with dense Cholesky, runs PCG, and independently checks ||Ax-b||/||b||. Compare its iteration count and solve time with the scalar-Jacobi result before changing the Auto escalation policy.

Two-level aggregation diagnostic

If block Jacobi improves the residual but still leaves a large iteration count, use the development two-level benchmark. It combines the SPD block-Jacobi term with a geometry-free Galerkin coarse correction,

M^-1 = B^-1 + Z (Z^T A Z)^-1 Z^T.

Z has one piecewise-constant mode per displacement component and contiguous node aggregate. For the RCM-ordered L-angle export, start with 1024 nodes per aggregate:

.\\benchmarks\\scripts\\bench-fem-twolevel.ps1 D:\\path\\to\\L-angle-K.mtx `
  -DofsPerNode 3 `
  -AggregateNodes 1024 `
  -Tolerance 1e-8 `
  -MaxIterations 3000

This is still a diagnostic path rather than the automatic HyBIT policy. Its purpose is to determine whether the remaining slow convergence is dominated by global smooth modes. Smaller aggregates create a larger, more expensive coarse space; larger aggregates create a cheaper but weaker coarse correction.

Rigid-body coarse-space experiment

For 3-D structural matrices with three displacement DOFs per retained node, bench-fem-rigid.ps1 uses six rigid-body modes per aggregate: Tx, Ty, Tz, Rx, Ry, Rz. The companion .coords file must be generated by the matching mf_solver HyBIT exporter so its RCM/reduced node order exactly matches the Matrix Market rows.

.\benchmarks\scripts\bench-fem-rigid.ps1 D:\path\L-angle-K.mtx -AggregateNodes 1024 -Tolerance 1e-8 -MaxIterations 3000

Using 1024 nodes/aggregate gives approximately the same coarse dimension as the 512-node translation-only experiment, isolating the value of elasticity-aware rotation modes from the effect of merely increasing coarse dimension.

Physical-load rigid-body benchmark

Once the matching reduced RCM-order load vector has been exported by mf_solver, pass it to the same rigid-body benchmark with -Rhs:

.\benchmarks\scripts\bench-fem-rigid.ps1 D:\path\L-angle-K.mtx `
  -Coordinates D:\path\L-angle-K.coords `
  -Rhs D:\path\L-angle-b.txt `
  -AggregateNodes 1024 `
  -Tolerance 1e-8 `
  -MaxIterations 3000

When a physical RHS is supplied there is no known x_exact=1, so the benchmark reports the solver and independently recomputed residual but omits the synthetic relative solution-error line.

Structural Auto API (0.6-r8)

For 3-D structural systems with node coordinates, HybitSolver now provides an explicit geometry-aware path. StructuralOptions chooses an aggregate size from a target dense coarse dimension, while solve_structural_csr32 and prepare_structural_csr32 apply the six-mode rigid-body two-level preconditioner. The normal generic solve_csr32 API remains unchanged.

Example:

.\benchmarks\scripts\bench-fem-structural-auto.ps1 L-angle-K.mtx `
  -Coordinates L-angle-K.coords `
  -Rhs L-angle-b.txt `
  -TargetCoarseDimension 1536 `
  -Tolerance 1e-8 `
  -MaxIterations 3000

For the 119355-node L-angle reduced system, the default target selects 512 nodes per aggregate and a 1404-dimensional rigid-body coarse space.

Graph-connected structural aggregation (0.6-r10)

The r9 baseline uses fixed-width aggregates in reduced RCM node order. r10 adds an experimental graph mode that infers node adjacency from nonzero 3x3 structural blocks and grows deterministic breadth-first connected aggregates. Use the same target coarse dimension for a controlled comparison:

.\benchmarks\scripts\bench-fem-structural-auto.ps1 L-angle-K.mtx `
  -Coordinates L-angle-K.coords `
  -Rhs L-angle-b.txt `
  -TargetCoarseDimension 1536 `
  -Aggregation graph `
  -Tolerance 1e-8 `
  -MaxIterations 3000

For a direct A/B run, use bench-fem-structural-aggregation.ps1; it executes the contiguous and graph policies sequentially with identical matrix, RHS, tolerance, and target coarse dimension. Keep the contiguous result as the reference. Prefer graph aggregation only if it reduces iterations or total solve time enough to justify its additional preparation cost.

L-angle Structural Auto development result (0.6-r11/r12)

For the exported 358065-DOF L-angle reduced stiffness matrix with the physical load RHS, target_coarse_dimension=1536 produced the following development measurements:

• contiguous rigid-body aggregation: 234 aggregates, coarse dimension 1404, 1370 PCG iterations;
• graph-connected rigid-body aggregation: 233 aggregates, coarse dimension 1398, 220 PCG iterations.

Both runs converged below 1e-8 verified relative residual. These figures are specific to this matrix, RHS, machine, and implementation revision and are not general speed claims.

Packed structural coarse factor experiment (0.6-r17)

The validated Graph/Additive structural path keeps the same Galerkin coarse operator and PCG algorithm, but stores its dense Cholesky lower factor in row-wise packed triangular form. A coarse dimension m therefore retains m(m+1)/2 factor values instead of m^2. Row offsets are precomputed so the packed indexing does not add integer address arithmetic to the inner forward and backward triangular-solve loops.

For the L-angle Graph sweep, the expected coarse-factor storage changes are approximately 14.91 MiB -> 7.46 MiB at actual coarse dimension 1398 and 59.39 MiB -> 29.70 MiB at actual coarse dimension 2790. Re-run the 1536 and 3072 target points before drawing a wall-time conclusion; the mathematical preconditioner is unchanged, so the iteration counts should remain identical up to floating-point roundoff.

Structural kernel profile and parallel CSR result (0.6-r18/r19)

The L-angle Graph/Additive profile at target coarse dimension 1536 measured about 24.39 ms per PCG iteration. The standalone kernel estimates attributed roughly 17.25 ms/call to the 28.24M-nnz CSR SpMV and 4.86 ms/apply to the complete rigid-body two-level preconditioner. This identified fine-grid SpMV as the dominant CPU kernel.

The r19 thread sweep then compared the validated serial CSR kernel with Rayon row parallelism without changing the preconditioner or convergence target. On the Ryzen 7 7800X3D development machine, 4 Rayon threads was the best tested full-PCG point: serial SpMV was about 13.24 ms/call, parallel SpMV about 5.09 ms/call, and the 220-iteration PCG solve improved from about 4.44 s to 2.73 s. Eight and sixteen threads gave similar SpMV timing but slightly worse full-PCG time, consistent with a memory-bandwidth-limited sparse kernel. These measurements are machine-specific.

Integrated structural SpMV policy (0.6-r20)

The public Structural Auto/prepared path now accepts StructuralSpmvPolicy. Auto keeps small problems serial and enables the validated parallel CSR operator for CSR32 systems with at least 1,000,000 nonzeros. The generic non-structural solve_csr32 path remains unchanged. For reproducible CPU tuning, the PowerShell structural wrappers accept -RayonThreads; the library does not hard-code the 4-thread optimum observed on the development machine.

Use the end-to-end A/B wrapper to verify the integrated path:

.\benchmarks\scripts\bench-fem-structural-auto-spmv.ps1 L-angle-K.mtx `
  -Coordinates L-angle-K.coords `
  -Rhs L-angle-b.txt `
  -TargetCoarseDimension 1536 `
  -Aggregation auto `
  -RayonThreads 4 `
  -Tolerance 1e-8 `
  -MaxIterations 3000

Parallel rigid-body preconditioner experiment (0.6-r21)

The integrated r20 L-angle benchmark preserved 220 iterations while reducing solve time from about 4.45 s with serial CSR to about 2.73 s with StructuralSpmvPolicy::Auto selecting parallel CSR at four Rayon threads. The next remaining CPU target is the rigid-body preconditioner itself.

The r21 experiment leaves the production Structural Auto policy unchanged and compares two preconditioners under the same parallel CSR operator. The experimental view executes independent 3x3 block-Jacobi factors, aggregate-local restriction, and node-local prolongation in parallel. The packed coarse triangular solve remains serial. An aggregate-to-node index is prepared once so restriction requires no atomics or locks.

.\benchmarks\scripts\bench-fem-structural-precond-parallel.ps1 L-angle-K.mtx `
  -Coordinates L-angle-K.coords `
  -Rhs L-angle-b.txt `
  -Threads 2,4,8,16 `
  -TargetCoarseDimension 1536 `
  -Tolerance 1e-8 `
  -MaxIterations 3000

Use full PCG wall time, not the isolated preconditioner speedup, to decide whether this should be integrated. Both PCG runs use the same parallel CSR SpMV so the comparison isolates preconditioner execution.

Structural parallel execution (0.6-r22 development)

Structural Auto can independently resolve fine-grid SpMV and rigid-body preconditioner execution policies. On large CSR structural systems, Auto selects parallel CSR SpMV and parallel block-Jacobi/restriction/prolongation kernels; the packed dense coarse triangular solve remains serial. The aggregate-to-node index used by the parallel restriction is cached during prepare and reused by solve-many.

Integrated PCG vector policy (0.6-r25)

The r24 L-angle sweep showed that parallel/fused PCG vector kernels are beneficial once the shared Rayon pool has enough workers: at 8 workers the 220-iteration experimental solve changed from about 1.998 s to 1.790 s, while the 2-worker path was slower. Structural Auto r25 therefore requires both a sufficiently large vector (n >= 131072) and at least four Rayon workers before enabling parallel/fused vector kernels automatically.

The integrated r25 production-path A/B at 8 workers held Graph aggregation, parallel CSR SpMV, and the parallel rigid-body preconditioner fixed. Serial PCG vectors took 1.900 s; Auto selected Parallel PCG vectors and took 1.726 s. Both paths converged in 220 iterations. The Auto run independently verified ||Ax-b||/||b|| = 9.378557e-9; total analysis+prepare+solve was 2.797 s. These are machine-specific development measurements.

Use bench-fem-structural-auto-pcg.ps1 to reproduce the Serial versus Auto comparison. The r25 structural production path is feature-frozen after this checkpoint; the next 0.6 work is release-gate validation.

HyBIT development status

Last updated: 2026-10-07

Published release

• Current published release: HyBIT 0.8.0.
• Release tag: v0.8.0.
• Six Rust crates are published as 0.8.0 on crates.io.
• hybit-ffi remains repository-only.
• Rust 1.73 remains the declared MSRV.

The 0.8 release freezes the validated GeneralSquare FGMRES/ILU(0), ABTM G1-G7, prepared execution, Rayon, and C/C++/Fortran thread-interoperability work.

0.8.0 release checkpoint

Current validated development sequence is complete through G7. The F2 base commit is 782ac38; F3 adds prepared multi-RHS Natural/RCM ILU(0) reuse measurement and documentation on top of that checkpoint.

Execution architecture

Completed CPU checkpoints include common problem/execution/backend policy, resident Krylov abstraction, serial resident PCG, Rayon vector execution, parallel CSR SpMV, parallel resident Jacobi, and structural resident cross-check.

GPU is an architectural target, not a production target here.

GeneralSquare E1-F7

1. E1 — restarted FGMRES with reusable V/Z workspace.
2. E2 — prepared GeneralSquare routing.
3. E4a — staged restart escalation.
4. E4b — budget-aware restart controller.
5. E5 — prepared canonical ILU(0).
6. F1 — Natural/RCM ILU(0) ordering sensitivity.
7. F2 — real nonsymmetric corpus preflight and larger-case evidence.
8. F3 — repeated prepared multi-RHS ordering/reuse and end-to-end amortization.
9. F4 — ILU(0) application/topology profiling; per-level Rayon triangular apply rejected.
10. F5 — explicit unsuitable-ILU fallback with retained Identity safety path and real missing-diagonal validation.
11. F6 — Natural/RCM ordering-selection signals, paired short probes, and amortized policy replay; no automatic production promotion.
12. F7 — Jacobi -> ILU(0) promotion study; strong ILU wins and strong regressions observed, with no validated automatic production selector.

Jacobi + fixed restart 30 remains default.

E5 ILU(0) canonicalizes unsorted/duplicate CSR, introduces no fill, stores diagonal positions as u32, uses a selective 1e-12 row-relative factor-pivot floor, and is reused across solve-many.

The public hard-B cross-check reproduced 276 iterations, true relative residual 9.455610e-9, 4.238 MiB ILU state, 4.500 MiB FGMRES(3) workspace, 8.738 MiB total persistent state, and zero adjusted pivots.

F1-F4 ordering/kernel work now shows three distinct real-nonsymmetric regimes: sherman5 with no measured RCM convergence benefit, raefsky3 with robust cross-RHS RCM benefit, and venkat25 with RHS-dependent convergence changes and only a small repeated 5-RHS end-to-end timing difference. Bandwidth and a single RHS are therefore insufficient automatic-reordering signals.

These are regression measurements, not universal performance claims.

Current routing limitations

• Spd -> PCG/Hybrid/structural.
• GeneralSquare -> FGMRES on Auto/Cpu.
• SymmetricIndefinite -> recognized; MINRES not implemented.
• no rectangular LSQR/LSMR;
• no complex arithmetic;
• ILU(0) triangular solve is serial;
• C ABI reports FGMRES/ILU0 codes but does not yet expose GeneralSquare configuration setters.

Published 0.7 validation

The exact 0.7 release source passed source-integrity, metadata, formatting, Clippy, workspace release tests, Rust 1.73 MSRV, C/C++/Fortran ABI/runtime, package/dry-run, and physical L-angle gates.

The L-angle release regression used 358065 free DOFs and 28239653 CSR nonzeros, a 1398-dimensional rigid-body coarse space, 220 iterations, and independently verified relative residual 9.378557e-9.

Next validation focus

GeneralSquare robustness is now closed through F7 for the current 0.8 development checkpoint.

Production semantics remain explicit:

• Jacobi is the GeneralSquare default;
• canonical ILU(0) remains opt-in;
• Ilu0Fallback remains the explicit missing-diagonal safety policy;
• Natural/RCM ordering remains explicit rather than automatically selected.

ABTM G1-G7 are validated through topology algebra, metadata-first pruning, region/local extraction, explicit ABTM ILU(0), fixed-size 3x3/6x6 block execution, explicit prepared restricted/local CSR operators, and explicit prepared Rayon execution. G5 retains an explicit DenseBlockCsrOperator path without automatic CSR/B3/B6 routing. G6 retains ordinary CSR for scalar SpMV and materializes AM / RA*R as compact prepared CSR. G7 exposes full-Rayon and task-limited prepared execution, keeps global-to-local gather serial, and promotes no hardware-specific automatic nnz selector. The next checkpoint is G8 distributed topology / partition / halo and MPI-host integration. GeneralSquare automatic-selection research can be revisited later with a broader held-out corpus or application-provided solve-horizon/cost information.

Do not add ILUT, fill, pivoting, or hidden automatic preconditioner changes as part of the F7 conclusion.

ABTM G1

G1 validates a numerical-value-independent sparse-of-64-bitmaps topology layer with AND/OR/AND-NOT/XOR, popcount, rank/select, invariants, occupancy telemetry, and general Boolean merge semantics.

The ten-matrix corpus shows that bitmap topology metadata is strongly compact on several matrices but not universal: low-occupancy cases can be break-even or worse than CSR metadata. Logical topology is therefore accepted while physical prepared layout selection remains adaptive and workload-specific.

Word-local rank/select is the intended packed-value addressing primitive. Row-wide rank/select remains a convenience/correctness API rather than a hot numeric-kernel path.

G2-G7 are now complete. Next: G8 distributed topology / MPI-host integration.

ABTM G2a

G2a adds a scalar reference prepared value stream aligned with the G1 logical topology and a sparse-row-dot kernel that intersects topology with a DofMask before loading numerical matrix values.

The first benchmark varies active support from 100% to 10% and records exact candidate/executed/skipped products, empty-word ratio, numerical agreement, and scalar timing versus unpruned CSR. This is an evidence checkpoint only; no production solver/backend routing changes.

ABTM G2c

G2c adds an experimental dual structural topology with row and column orientations but no duplicated numerical values. The benchmark compares bitmap-word row/column support intersection with an explicit sorted-index baseline on deterministic sparse-dot pairs.

This checkpoint is intended to decide the metadata representation for later matrix-product, ILU symbolic-intersection, and coarse/Galerkin experiments.

ABTM G2d

G2d extends the dual-topology experiment to numerical row-by-column sparse dot products. It compares an explicit CSR/CSC-like baseline with two bitmap variants: one numerical copy plus a column source-index map, and duplicated column values.

The goal is to determine whether dual structural metadata can remain value-single-copy in practice or whether hot transpose-oriented kernels require a second packed numerical stream.

ABTM G2e

G2e evaluates numerical row/column products with the existing adaptive Sparse/Bitmap/Dense tile classification in both orientations. Sparse tile pairs avoid per-product rank, dense tiles use direct offsets, and bitmap tiles retain mask/rank addressing.

The purpose is to determine whether the G2d regressions on low/intermediate occupancy are a fixed-bitmap execution artifact before changing thresholds or promoting a numerical backend.

ABTM G2f

G2e confirms that adaptive Sparse/Bitmap/Dense execution is the correct direction, but also exposes that the current fixed 16-byte tile descriptor is not a compact Sparse physical encoding.

G2f benchmarks a packed variable-payload representation in which Sparse tiles store byte offsets and Bitmap/Dense tiles store 64-bit masks. This is intended to test the physical metadata model before threshold sweeps or production backend promotion.

ABTM G2g

G2g replaces the G2f byte-oriented payload decoder with typed compact streams and an 8-byte descriptor while retaining adaptive Sparse/Bitmap/Dense execution. It tests whether most of G2f’s storage reduction can be retained without losing the G2e numerical performance characteristics.

ABTM G2 closed

G2 metadata-first product pruning is validated and closed. The final CPU candidate is the G2e-style adaptive Sparse/Bitmap/Dense numerical execution with dual row/column structural preparation and explicit fallback.

The ten-matrix G2e corpus contains both clear ABTM wins and clear explicit wins, so no universal replacement or production selector threshold is promoted. The complete evidence is in ABTM_G2_CLOSEOUT.md.

Next checkpoint: G3 region growth, overlap/multiplicity, and local submatrix extraction.

ABTM G3a region growth

G3 is active. G3a adds deterministic k-hop region growth over support(A) union support(A^T) using the dual ABTM topology and an exact CSR-plus-transpose reference cross-check.

The checkpoint measures frontier nodes, topology words, candidate neighbor bits, explicit neighbor entries, preparation cost, and scalar execution time. Overlap/multiplicity and local submatrix extraction remain subsequent G3 work.

ABTM G3b overlap and multiplicity

G3b adds exact per-node multiplicity over G3a-grown regions, overlap-node statistics, pairwise overlap accounting, and an end-to-end growth-plus-multiplicity benchmark. Region equality is still verified against CSR + transpose-CSR before overlap metrics are accepted.

G3 local submatrix extraction remains the next step after the overlap corpus is validated.

ABTM G3c structural local-submatrix extraction

G3c adds exact structural extraction of A[R,R] from a region mask and dual ABTM topology. Extracted global-node lists, row pointers, and local column indices are cross-checked against direct CSR extraction for every benchmark region.

Numerical value gathering is intentionally deferred until structural pruning, mapping scratch cost, and extraction performance are measured.

ABTM G3d prepared local numeric refresh

G3d evaluates symbolic/numerical separation for local submatrices. A prepared plan binds the G3c local structural pattern to source CSR value positions and supports repeated numerical refresh without rebuilding the region or local pattern.

The benchmark reports one-time plan preparation, refresh-only time, direct full-extraction time, plan memory, numerical agreement under changed values, and the measured refresh count needed to amortize preparation.

ABTM G3 closeout

G3 region operations are validated and closed for this 0.8 checkpoint.

• G3a: deterministic undirected structural region growth;
• G3b: overlap and multiplicity with pair-overlap identity checks;
• G3c: exact structural A[R,R] extraction;
• G3d: prepared local numerical refresh under fixed topology.

The primary architectural result is a separation of concerns: ABTM is strongest as a reusable symbolic/topology layer, while materialized local sparse/numeric storage should remain evidence-driven. One-shot local pattern extraction did not show a universal speed advantage, whereas repeated prepared value refresh amortized its plan in roughly 4–8 refreshes on the six-matrix corpus when the dual topology was already available.

nd3k remains a deliberate high-overlap memory stress case and should not drive an unbounded all-regions-at-once policy.

Next work is G4 ABTM symbolic/numeric ILU(0) comparison against canonical CSR ILU(0).

ABTM G4a ILU(0) symbolic intersection

G4 is active. G4a compares canonical CSR ILU(0) symbolic target lookup against prepared ABTM row-word intersection while preserving identical canonical structure.

The first gate requires exact per-pivot target-list agreement. Numerical factorization, pivot handling, triangular apply, and solver residual behavior remain unchanged until symbolic evidence is established.

ABTM G4b ILU(0) numeric intersection

G4b evaluates numeric ILU(0) updates driven directly by ABTM word intersections. A per-word structural prefix maps each intersection bit to the existing canonical CSR factor position, retaining conventional numerical storage and avoiding per-product plan memory.

The gate requires factor and triangular-apply agreement against the canonical CSR reference before timing. Production routing remains unchanged.

ABTM G4c ILU(0) rank-LUT addressing

G4c isolates the G4b numerical-addressing overhead by replacing per-update rank-by-popcount with an all-word u8[64] rank table. CSR, G4b popcount, and G4c LUT factors are cross-validated before timing.

The experiment reports LUT preparation time and bytes separately. No production selection rule is introduced at this stage.

ABTM G4d adaptive rank-LUT sweep

G4d follows the successful G4c all-word LUT experiment with an occupancy-threshold sweep. The goal is to determine whether most of the numerical speedup can be retained with substantially less LUT memory.

Every threshold is cross-checked against canonical CSR factors and triangular apply before timing. Production ILU(0) remains unchanged.

ABTM G4f explicit production ILU(0)

G4f promotes the validated direct rank-LUT factorization algorithm into an explicit production constructor and GeneralSquare policy. The factor object retains only the canonical CSR factor arrays; ABTM topology and rank LUTs are temporary preparation metadata.

Automatic matrix-level routing is intentionally deferred. G4e held-out results are encouraging, but selector cost and a cheap production predictor still need separate validation.

ABTM G4 closeout

G4 is closed. The explicit production Ilu0Abtm path preserves canonical CSR factor storage and triangular application while using temporary ABTM topology and direct word-rank LUTs during factor construction.

The eleven-matrix production constructor corpus was exact on every supported case and split 6 ABTM wins / 5 CSR wins. Strong wins include nd3k (~5.28x), x104 (~2.34x), cant (~2.33x), and raefsky3 (~2.06x).

The production corpus also showed that the earlier permissive symbolic selector must not be promoted: venkat25 and cfd1 have favorable symbolic compression but lose after full ABTM preparation cost. Automatic matrix-level routing is therefore deferred. G5 block-ABTM is next.

ABTM G2-G5 closeout

G2 validates metadata-first support intersection/pruning and establishes that logical row/column topology can drive candidate elimination before numerical value access.

G3 validates topology-driven region growth, overlap/multiplicity, local structural extraction, and repeated local numeric refresh. The architectural boundary remains deliberate: ABTM discovers/prunes reusable structure while conventional prepared local numerical layouts remain appropriate downstream.

G4 validates ABTM symbolic ILU(0) intersection and exposes explicit Ilu0Preconditioner::from_csr32_general_abtm / GeneralSquarePreconditionerPolicy::Ilu0Abtm. The returned factor remains the canonical CSR ILU(0) representation; ABTM preparation metadata is scratch. Automatic ILU route selection remains deferred.

G5 validates node-oriented dense block execution only where structure supports it. Generic block traversal was rejected. Fixed-size 3x3/6x6 tail-specialized kernels produced substantial CSR speedups on high-fill development cases, while low-fill matrices remained better on CSR. A frozen held-out selector classified five matrices correctly with zero false-positive block routes; inline_1 provided a positive B3 held-out case at about 1.78x CSR kernel speed and about 0.70x CSR storage.

The production G5 result is therefore explicit rather than automatic: DenseBlockCsrOperator with DenseBlockSize::B3 / B6 is public and implements LinearOperator, but existing scalar CSR and solver defaults are unchanged. Positive held-out evidence for B6 is still missing, so automatic CSR/B3/B6 routing is not promoted.

Next: G6 ordinary and masked/restricted SpMV.

HyBIT roadmap

This roadmap describes direction, not guaranteed dates.

Published baseline

HyBIT 0.7.0 is published at immutable tag v0.7.0, commit 1fdcd6a1b8127c84306c38c3fdbad42563538ad8.

Completed in 0.7

• generic algebraic two-level SPD/PCG correction;
• graph/contiguous aggregation and smoothed transfer;
• parallel coarse transfer and storage policies;
• coarse-first resumable PCG;
• prepared coarse reuse;
• validated structural Graph rigid-body parallel CPU path;
• release integrity/MSRV/package/ABI/real-FEM gates.

Completed on develop/0.8.0 through G7

Execution architecture

• problem/execution/backend/solver/preconditioner policy architecture;
• resident Krylov abstraction;
• serial/Rayon resident PCG validation;
• parallel CSR/Jacobi execution;
• structural resident cross-check.

GeneralSquare

• restarted FGMRES with separate V and Z;
• reusable FGMRES workspace;
• prepared GeneralSquare routing;
• fixed, escalating, and budget-aware restart;
• canonical prepared ILU(0);
• selective row-relative pivot stabilization;
• prepared ILU reuse and memory telemetry;
• F1 Natural-versus-RCM ordering harness with original-system residual verification;
• real nonsymmetric ordering evidence showing both a strong RCM benefit (raefsky3) and a no-benefit case (sherman5), ruling out unconditional RCM promotion;
• F2 multi-matrix preflight/corpus harness with machine-readable CSV evidence;
• larger real nonsymmetric venkat25 evidence showing a moderate RCM ILU(0) benefit while Jacobi remains ordering-invariant;
• explicit corpus evidence that structurally missing diagonals are a common current-applicability boundary;
• F3 prepared Natural/RCM ILU(0) multi-RHS harness with solver-only and end-to-end permutation accounting;
• repeated 5-RHS evidence showing robust cross-RHS RCM benefit on raefsky3 but only small/timing-sensitive aggregate benefit on venkat25;
• F4 serial triangular-apply/dependency profiling, with per-level and width-threshold Rayon prototypes rejected because synchronization/scheduling overhead outweighed available level parallelism.

Immediate next work

GeneralSquare robustness

GeneralSquare robustness work is closed through F7 for this 0.8 checkpoint.

1. Keep Jacobi as the default GeneralSquare preconditioner.
2. Keep canonical ILU(0) explicit/opt-in.
3. Keep Ilu0Fallback explicit for the missing-diagonal safety path.
4. Keep Natural/RCM ordering explicit; F6/F7 evidence does not justify a universal automatic selector.
5. Revisit automatic selection only with broader held-out validation or application-provided solve horizon / cost budget.
6. Treat drop-tolerance/fill, pivoting, and multilevel ILU as separate future preconditioners rather than silently changing canonical ILU(0).

ABTM topology and metadata-first execution

After the current GeneralSquare robustness sequence, reintroduce ABTM as a symbolic/topology layer rather than assuming it must replace CSR everywhere.

1. G1 (validated): bitmap topology algebra: AND/OR/AND-NOT/XOR, popcount, rank/select, invariants.
2. G2 (validated): metadata-first product pruning with candidate/executed work metrics.
3. G3 (validated): region growth, overlap/multiplicity, local extraction, and symbolic/numeric reuse boundaries.
4. G4 (validated): ABTM ILU(0) symbolic intersection and explicit production ABTM factorization path; automatic routing deferred.
5. G5 (validated): 3x3/6x6 block-topology characterization and explicit fixed-size dense block-CSR numerical operator; automatic CSR/B3/B6 routing deferred pending broader positive held-out evidence.
6. G6 (validated): ordinary scalar ABTM SpMV remains diagnostic; fixed column restrictions and graph-local restrictions use explicit prepared compact CSR, with direct or ABTM-assisted preparation selected explicitly by workload.
7. G7 (validated): explicit full-Rayon and task-limited prepared CSR execution; serial gather retained; automatic hardware-specific size routing deferred.
8. G8 (next): distributed topology, partition/halo extraction, and MPI-host integration.
9. G9: GPU/CubeCL-specific prepared execution after the CubeCL API stabilizes.

Track occupancy, metadata bytes, rank/select cost, bitmap operations, candidate/executed products, pruning ratio, preparation cost, and execution time. Backend promotion remains evidence-driven.

Documentation/API completeness

• maintain separate user, numerical, architecture, and development docs;
• keep examples runnable under gates;
• add GeneralSquare selection to C/C++/Fortran only after the Rust policy is stable enough;
• document report fields by solver path.

Subsequent numerical work

• MINRES for SymmetricIndefinite;
• genuine local spectral/GenEO-style SPD enrichment;
• symbolic/topological reuse separated from numeric refactorization, with the ABTM topology track providing the structural representation experiments;
• rectangular LSQR/LSMR later;
• complex arithmetic later.

Parallel and large-scale work

• parallel setup where worthwhile;
• improve CPU sparse scheduling and NUMA/cache behavior;
• GPU-resident execution after the selected CubeCL interface is stable;
• GPU-suitable preconditioners rather than blindly porting serial ILU solves;
• distributed-memory domain decomposition;
• out-of-core state only when problem scale justifies it.

Release discipline

Keep release tags immutable, preserve the C ABI where practical, update language bindings together, and retain source-integrity, formatting, Clippy, MSRV, package, ABI/language, and numerical regression gates.

Non-goal

HyBIT does not claim universal sparse-solver coverage or universal performance improvement. Automatic policy changes require evidence across representative matrices.

Publishing HyBIT

Public repository: https://github.com/michioga/hybit

User-facing Rust crate: https://crates.io/crates/hybit

HyBIT 0.8.0 release procedure

HyBIT 0.8.0 publishes the same six Rust crates as 0.7.0. hybit-ffi remains repository-only but its C ABI, C++ wrapper, Fortran module, and OpenMP/Rayon thread-interoperability path are release-qualified.

The 0.8 release candidate must pass:

• source integrity;
• exact 0.8.0 workspace/dependency metadata;
• formatting and Clippy;
• workspace release tests and all targets;
• Rust 1.73 MSRV;
• C/C++/Fortran ABI/runtime;
• OpenMP/Rayon environment/API interoperability;
• crates.io package checks and dependency-ordered dry-runs;
• the physical L-angle structural regression and prepared reuse;
• mdBook build.

After qualifying the exact source commit on develop/0.8.0, move that exact source to main and rerun the complete gate there. Publish crates in dependency order, create immutable tag v0.8.0, and create the GitHub Release from RELEASE_NOTES_0.8.0.md.

GitHub Pages deploys the frozen 0.8 documentation from main.

Threading integration and oversubscription requirements are documented in THREADING.md.

HyBIT 0.7.0 publication record

HyBIT 0.7.0 was published on 2026-09-27 from the exact validated source commit:

1fdcd6a1b8127c84306c38c3fdbad42563538ad8

The immutable release tag is v0.7.0, and the GitHub Release uses RELEASE_NOTES_0.7.0.md.

The published Rust crates are:

1. hybit-core
2. hybit-matrix
3. hybit-krylov
4. hybit-precond
5. hybit-auto
6. hybit

hybit-ffi remains repository-only and is not published to crates.io. Its C ABI, C++ wrapper, and Fortran ISO_C_BINDING consumer path are built and runtime-tested by the repository release gate.

0.7.0 qualification record

The r32 numerical production line was frozen before release. Later watchdog, local-energy-gate, and filtered spectral-enrichment experiments were excluded from 0.7.0.

The exact release commit passed the complete release-candidate-gate.ps1 on develop/0.7.0, was fast-forwarded to main, and passed the complete gate again on main before publication. In the immutable v0.7.0 tree that script is at repository root; the post-release 0.8 development tree retains the 0.7 qualification scripts under tools/release/0.7/.

The qualifying gate covered:

• source-integrity verification, including tracked-file/MANIFEST.txt agreement;
• workspace metadata and version checks;
• formatting and Clippy with warnings denied;
• workspace release tests and all release targets;
• Rust 1.73 MSRV;
• Rust/C ABI plus C, C++, and Fortran build/runtime checks;
• package inspection and crates.io dry-run validation;
• the real L-angle structural regression;
• independently verified residual, iteration guard, and prepared solve-many reuse.

The release-reference L-angle solve converged in 220 iterations with independently verified relative residual 9.378557e-9. Prepared reuse also passed.

Immutable release rule

Published crate versions and release tags are immutable.

Do not move v0.7.0 to a later documentation or development commit. A serious defect in 0.7.0 must be corrected with a new version such as 0.7.1 or a later release, not by rewriting the existing tag or crates.io artifacts.

Post-release documentation commits may advance main; they do not alter the 0.7.0 release source.

crates.io authentication

Authenticate locally with a crates.io API token using cargo login. Treat the token as a secret; never commit Cargo credentials or paste tokens into logs.

For GitHub Release automation, GitHub CLI may be authenticated with gh auth login. Authentication tokens and device codes must not be committed or included in release logs.

Dependency order for future workspace releases

HyBIT is a Cargo workspace. Higher-level crates depend on lower-level crates, so publication must preserve this dependency order:

hybit-core
    |
    +--> hybit-matrix
    +--> hybit-krylov
              |
              v
         hybit-precond
              |
              v
          hybit-auto
              |
              v
            hybit

hybit-matrix and hybit-krylov may be published at the same dependency level after hybit-core is visible in the crates.io index.

For each future release:

1. Freeze the exact numerical source set and exclude unintended experimental/generated artifacts.
2. Regenerate MANIFEST.txt and SOURCE_SHA256.txt.
3. Run the complete release-candidate gate without release-qualifying skip switches.
4. Merge the exact validated commit to main.
5. Rerun the complete gate on main.
6. Publish crates in dependency order, confirming each dependency level is visible in crates.io before publishing dependents.
7. Create an immutable annotated version tag on the exact published source commit.
8. Create the GitHub Release from the matching release notes.
9. Perform any post-release documentation synchronization in a later docs-only commit without moving the release tag.

The release gate records performance measurements for regression context, but correctness, residual verification, compatibility, and prepared reuse remain the publication criteria.

ABTM G3 closeout: regions, overlap, local extraction, and symbolic reuse

G3 validates ABTM as a structural/topology layer for region-oriented sparse operations. The checkpoint covers deterministic region growth, overlap and multiplicity, structural extraction of A[R,R], and prepared numerical refresh when topology remains fixed.

Status

G3 is validated and closed for the 0.8 development checkpoint.

The results do not justify replacing CSR universally. They support a split architecture:

• use ABTM dual topology for structural graph operations and metadata-first pruning;
• materialize ordinary local sparse structure when downstream numerical kernels require it;
• cache symbolic/numerical address plans when the same local structure is refreshed repeatedly.

G3a: region growth

G3a grows seed regions over

support(A) union support(A^T)

using a frontier formulation. Every ABTM region was cross-checked against an explicit CSR + transpose-CSR reference.

Two-hop corpus:

The one/two/three-hop sweep showed that tiny frontiers are dominated by fixed bitset and allocation costs. As frontier work grows, matrices with many structural entries represented per topology word benefit strongly. nd3k became a saturation stress case: at three hops the average region was about 4612 nodes out of a 9000-node universe.

This confirms that topology-word compression can translate directly into graph traversal savings, but only when enough work exists to amortize setup and frontier bookkeeping.

G3b: overlap and multiplicity

For grown regions R_i, G3b defines

m(v) = number of regions containing node v.

Correctness was checked three ways:

• every ABTM-grown region equals the CSR reference region;
• every per-node multiplicity equals a direct reference count;
• pair overlap obeys sum_{i<j} |R_i intersect R_j| = sum_v C(m(v), 2).

Two-hop corpus:

Multiplicity reduction itself was cheap compared with region growth, so the pipeline largely inherited G3a behavior. nd3k is intentionally treated as a high-overlap stress case: all 9000 nodes were covered and overlapped, 1594 of 2016 region pairs intersected, and maximum multiplicity reached 31.

G3c: structural local-submatrix extraction

G3c materializes the structural CSR-like pattern of A[R,R]. ABTM intersects row topology words with the region mask before enumerating retained columns. The result was checked exactly against direct CSR extraction:

global_nodes
row_ptr
col_idx

all matched for every tested region.

Extraction-time ratios:

The main conclusion is negative but useful: after the local structure must be fully materialized, output enumeration and global-to-local mapping dominate enough that topology-word compression alone is not a general speed advantage.

nd3k produced 38,158,139 local structural entries across the 64 regions, so the extraction benchmark was dominated by output construction despite a topology-word / CSR-entry ratio near 0.05.

G3d: prepared local numerical refresh

G3d tests the reuse case suggested by G3c. A plan is prepared once that binds each local structural entry to source global CSR value position(s). Later value refreshes reuse the local topology and address map without rebuilding regions, local numbering, row pointers, or column indices.

All six corpus matrices had:

mismatched_regions=0
mismatched_patterns=0
mismatched_values=0
max_scaled_error_original=0
max_scaled_error_refresh=0

for the tested values. All 64 plans on each corpus matrix used the direct single-source path; the Matrix Market loader had already canonicalized duplicates. Duplicate-source summation is covered separately by unit tests.

The reported break-even count amortizes plan preparation only. If the dual topology does not already exist, its preparation cost must also be included. In the intended G3/G4 architecture the dual topology is a reusable global symbolic object, so both views are relevant.

The nd3k stress case also exposes the memory boundary: 43,338,226 local structural entries across 64 heavily overlapping regions require about 497 MiB for the current plans before local numerical values themselves are materialized. Production code should avoid retaining unnecessary identity source-pointer arrays for direct-source plans and should not build all high-overlap region plans simultaneously without a memory policy.

Interpretation boundary

G3d demonstrates the value of prepared symbolic reuse. It does not prove that ABTM has an intrinsic advantage over every possible prepared CSR implementation: a CSR-based implementation can also cache local source positions once the structure is known.

The ABTM-specific evidence is strongest in G3a/G3b, where dual topology and word-level set algebra reduce structural graph work. G3c/G3d then show the appropriate interface boundary: use ABTM to discover and prune structure, but allow conventional local sparse/numeric layouts downstream.

The G3 benchmark programs use independent deterministic seed streams. Therefore absolute region totals from G3a, G3b, G3c, and G3d are reproducible within each checkpoint but are not intended for direct stage-to-stage equality comparisons.

G3 architectural decision

Promote the following direction:

1. dual ABTM topology for graph/region operations;
2. deterministic frontier growth over support(A) union support(A^T);
3. explicit overlap/multiplicity metadata;
4. exact materialization of local sparse patterns when required;
5. prepared source-address maps for repeated numerical refresh;
6. evidence-driven memory policy for large or highly overlapping region sets.

Do not promote:

• universal ABTM replacement of CSR local matrices;
• one-shot local pattern extraction as a performance claim;
• all-region materialization without a memory budget;
• G3d refresh ratios as proof that bitmap value storage is superior.

Next: G4

G4 should compare ABTM symbolic/numeric ILU(0) work against canonical CSR ILU(0), while preserving the G2/G3 separation:

• symbolic candidate intersection and pruning from topology;
• conventional prepared numerical factors where that wins;
• exact factor/apply cross-checks and original-system residual verification;
• preparation, memory, factorization, and repeated-apply accounting;
• no automatic promotion from a single matrix or occupancy threshold.