Keyboard shortcuts

Press ← or → to navigate between chapters

Press S or / to search in the book

Press ? to show this help

Press Esc to hide this help

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.