perf(engine): sparse dynamic solvers — time history and harmonic - #180
Open
diegokingston wants to merge 22 commits into
Open
perf(engine): sparse dynamic solvers — time history and harmonic#180diegokingston wants to merge 22 commits into
diegokingston wants to merge 22 commits into
Conversation
…te degrees Replace the exact elimination game (per-node HashSet adjacency, explicit fill-edge insertion) with the canonical AMD scheme (Amestoy-Davis-Duff 1996, same structure as CSparse's cs_amd): eliminated variables become quotient-graph elements whose external adjacency is kept explicitly, degrees are approximate external degrees, absorbed elements are dropped during degree updates, and variables whose live adjacency collapses into the new clique are mass-eliminated without touching the pivot heap. Determinism is preserved: all scans follow array insertion order and the pivot heap is keyed on (degree, index), so ties resolve to the smallest index. Simplifications vs full AMD: no supervariable merging and no post-ordering (the caller builds the elimination tree afterwards). On a 256x256 5-point grid Laplacian (n=65536), amd_order goes from ~3.0s to ~75ms (40x) and produces slightly less fill than the old exact-degree ordering (l_nnz 2.02M vs 2.33M). Fill quality on ~140 patterns (grids, arrowheads, banded, random SPD) is pinned by a property test against the old implementation kept as a test-only reference. The modal-shell and buckling-3D goldens in sparse_mass.rs are recaptured: the new elimination order changes FP rounding in the K factorization (~1e-9 relative), and the new values are as close or closer to both the analytical Euler values and the dense generalized Lanczos path. Also expose symbolic_cholesky_with_perm to allow injecting an explicit permutation into the symbolic factorization (used by the fill-quality property test).
…rmal plate stress - numeric_cholesky: guard as !(diag > threshold) so a NaN pivot (where every <= comparison is false) is rejected instead of returning a NaN-filled factor reported as success. - shell_benchmark: the restrained/free thermal plate benchmarks only checked displacements and reactions, which pass even when the thermal strain is not subtracted in stress recovery. Assert the actual stress contract: clamped plate reports -E*alpha*dT/(1-nu), free plate ~0.
Replace the left-looking merge (O(nnz(L)^2): every column contributing to
a row re-scanned per row entry) with the standard two-step construction:
1. Etree of the permuted matrix directly from the graph of A (Liu's
algorithm with path compression, O(nnz*alpha(n))).
2. Column patterns by union over etree children:
struct(L[:,j]) = {j} u struct(A[:,j]) u (u_c struct(L[:,c]) \ {c}),
so each column of L is scanned exactly once, by its parent.
The pattern and the tree are identical to the previous implementation,
bit for bit; a cfg(test) copy of the old merge serves as reference and
asserts exact equality of l_col_ptr, l_row_idx and parent on assembled
stiffness matrices from the nave-industrial, tower, space-truss and
building-case1 fixtures, under both AMD and RCM orderings.
Measured (release, symbolic_cholesky total, AMD included):
nave-industrial 77ms -> 5.5ms, building-case1 7.8ms -> 0.97ms.
The symbolic phase now partitions columns into fundamental supernodes (parent[j]=j+1 and column counts differing by one) and precomputes, per supernode, the list of prior supernodes that update it. The numeric phase factors each supernode as a dense cache-resident panel: scatter A, gather updates from prior supernodes through the dense trapezoidal layout that supernodes induce inside the CSC storage (column s+d holds row set R[d..]), dense partial Cholesky on the panel, scatter back. The output format (CSC L + permutation) is unchanged, so solves and all callers are unaffected. This also moves the per-call rebuild of the nonzero-column lists out of numeric_cholesky: the update structure is pattern-derived and now lives in SymbolicCholesky, which P-Delta and the Lanczos paths reuse across factorizations. Correctness is pinned by a cfg(test) copy of the simplicial algorithm: l_values must match it to 1e-9 relative on synthetic SPD matrices (dense, banded, arrowhead) and on the assembled nave-industrial, tower, space-truss and building-case1 stiffness matrices, under AMD and RCM. The modal-shell golden in sparse_mass.rs is recaptured (summation order changed, ~1e-9..1e-8 relative on lambda, same convention as the AMD recapture) and the test gains the dense-path cross-check the frame test already had. Measured (release, this branch vs main): - 256x256 5-point grid (n=65536): symbolic 3.58s -> 81ms (44x), numeric 177ms -> 110ms, nnz_L 2.33M -> 2.02M (less fill). - frame-like grid of 6-DOF cliques (n=9600): symbolic 977ms -> 26ms (38x), numeric 60ms -> 40ms; 1228 supernodes, max width 390. - nave-industrial fixture: symbolic 46ms -> 2.7ms, numeric 2.5 -> 1.5ms.
A single rigid diaphragm used to pull the whole constrained path into dense cubic territory: build_constraint_transform allocated C as dense n_total x n_independent (~370 MB at 6.8k DOFs), FreeConstraintSystem kept a dense nf x p C_ff, and C'KC ran as two dense triple loops (~2*nf^3 flops) over a dense K_ff that solve_constrained_3d extracted from a dense n x n assembly. C is near-identity: one unit entry per independent DOF plus a few entries per dependent one. This change keeps it that way end to end: - SparseTransform (dual CSR/CSC) replaces the dense c_matrix; the chained master-of-master substitution keeps its multi-pass semantics but works on sparse rows with a generation-stamped accumulator. - FreeConstraintSystem stores c_ff sparse; reduce_matrix (dense in/out, used by the nonlinear and dynamic paths) is now O(nnz(C)*nf) instead of O(nf^3), and reduce_matrix_sparse streams K's triplets into a sparse C'KC at O(nnz(K)*deg^2) — near-linear for MPC transforms. - solve_constrained_3d is sparse end to end: assemble_sparse_3d with k_full (prescribed-DOF correction via sparse_cross_block_matvec, reactions via sym_mat_vec), sparse C'KC, sparse Cholesky with the dense LU fallback kept for the singular case. The small/large path split (dense below 64 DOFs) is gone — one diagnostics contract for all sizes. - modal 2D/3D constrained paths reduce K and M as sparse triple products and use the sparse Lanczos; no dense mass assembly, no to_dense_symmetric undo. 2D linear constrained keeps its dense k_ff (the 2D sparse assembly has no k_full yet) but already benefits from the sparse C and the O(nnz(C)*nf) reduction. Measured on a 5x5x10 frame (2.2k free DOFs) with one rigid diaphragm per floor, on top of the quotient-graph AMD branch: 22.3 s -> 31.9 ms (~700x), bit-identical displacements and reactions. Correctness: a unit test pins reduce_sparse against reduce_dense and a naive triple product (regression for an ordered-slot double-count caught during development); the full suite (7101 tests) passes, including the chained-constraint depth tests and the constrained shell-stress parity. One diagnostics test updated: small constrained 3D models now report SparseCholesky since the dense small-model path no longer exists. Stacked on perf/amd-quotient-graph (#176).
Time history (2D/3D): assemble K and M as CSC, reduce constraints as sparse triple products, Rayleigh damping as a0*M + a1*K over the union pattern, K_eff via CscMatrix::linear_combination, factor once with sparse Cholesky. Dense LU fallback stays for small models, size-capped (MAX_DENSE_FALLBACK_DOFS) so a large model fails loudly instead of exhausting the WASM address space. Peak reactions keep a sparse full-K via assemble_sparse_2d_ex(build_k_full). Harmonic (2D/3D): the sparse modal path now covers constrained models too (sparse C'KC/C'MC reduction + reduced target DOF mapping) instead of always densifying; the dense modal/direct fallbacks and the 2n x 2n complex LU get the same dense-memory ceiling. New integration tests: harmonic with EqualDOF/RigidLink constraints on 2D (sparse branch) and 3D, including slave==master response parity and the dependent-target error.
The buckling paths built K_g as a dense n×n matrix and then converted it to CSC with from_dense_symmetric for the sparse Lanczos eigensolver — the n² allocation the sparse path exists to avoid. The constraint paths additionally densified K_ff and ran dense Jacobi. - geometric_stiffness: build_kg_from_forces_2d/3d and the five shell add_*_geometric_stiffness_3d now scatter into a KgTriplets sink that keeps only the free×free block and converts to lower-triangle CSC (duplicates summed). Element loops emit each unordered local pair once (lower triangle) — emitting the full symmetric square would double-count off-diagonals through from_triplets' duplicate summation (caught by the Euler buckling goldens during development). - buckling 2D/3D: no dense K_g, no extract/negate roundtrip; the constraint path reduces K and -Kg with the sparse triple product and uses the same sparse shift-invert Lanczos as the unconstrained path (dense Jacobi gone from buckling). Full suite green (7101), including the buckling goldens and the sparse-vs-Jacobi parity test (whose manual reference path now densifies the CSC triplets instead of the other way around). Stacked on perf/sparse-constraint-transform (#177).
numeric_cholesky used to call permute_symmetric on every factorization: tripletize + sort_unstable O(nnz log nnz) plus six nnz-sized temporaries, all of it purely symbolic work. The symbolic phase now stores the permuted structure (pa_col_ptr/pa_row_idx) and pa_src[p] = index into the original values array, so a numeric factorization permutes values with a single O(nnz) gather and no allocations beyond the output. The pattern-reuse contract is unchanged (callers already had to pass the same-structure matrix); the doc comment now states it explicitly. Measured on the 256x256 grid (n=65536, nnz_L=2.02M), 5 repeated numeric factorizations with the symbolic reused (the P-Delta pattern): ~116 ms -> ~100 ms per factorization (medians of 3 runs; noisy box). The win is bigger in allocation pressure than in wall time: no per-call triplet buffers. Stacked on perf/amd-quotient-graph (#176).
2D paths dispatch at SPARSE_THRESHOLD: small models keep the original dense code, large models assemble K (and M for Craig-Bampton) as CSC, reduce constraints as sparse triple products, keep K_II as a sparse principal submatrix, and factor it once with sparse Cholesky (dense LU fallback capped by MAX_DENSE_FALLBACK_DOFS). Craig-Bampton interior modes come from sparse shift-invert Lanczos. Guyan reactions use the sparse full-K matvec instead of dense K_rf/K_rr extraction. 3D paths stop densifying when constraints are present (reduce_matrix_sparse instead of reduce_matrix(to_dense)) and get the same sparse K_II factorization; Craig-Bampton 3D drops the dense O(n^2) mass assembly. Blocks involving only boundary DOFs stay dense (nb is small by design); the small-model 2D path is byte-identical to the old code.
Entries that are analytically zero pick up rounding-level asymmetry from the different operation order of the sparse products; comparing the difference against the entry pair itself fails on near-zero couplings. Scale against max|K| instead.
Above SPARSE_THRESHOLD the 2D/3D co-rotational Newton loops now convert the extracted free tangent block to CSC, reduce constraints with reduce_matrix_sparse, and factor with sparse Cholesky — the symbolic factorization is computed once per solve call and reused across iterations (pattern-fingerprinted; rebuilt only if the CSC structure changes). Modified NR caches the sparse factor per increment. Non-SPD tangents fall back to dense LU, capped by MAX_DENSE_FALLBACK_DOFS. Below the threshold the dense path is unchanged. New validation test: 30-element cantilever (nf = 90) converges through the sparse path, tip deflection within 1% of PL^3/3EI, dense/sparse and full/modified NR parity within 1e-4.
…ment control FactoredTangent gains a Sparse variant (numeric Cholesky over a cached symbolic factorization, fingerprinted by the CSC pattern and reused across all steps/iterations of a solve call). Above SPARSE_THRESHOLD the constraint reduction moves to reduce_matrix_sparse; non-SPD tangents fall back to dense LU capped by MAX_DENSE_FALLBACK_DOFS. Factor-once / solve-many is preserved (arc-length corrector still solves 2 RHS per factorization). Below the threshold the dense path is byte-identical. New validation tests: 30-element cantilever (nf=90, sparse) vs 21-element mesh of the same cantilever (nf=63, dense) — arc-length and displacement-control load factor and tip displacement agree within 2%.
…solvers SparseSymbolicCache / cached_symbolic / tangent_free_sparse / solve_tangent_sparse move from corotational.rs (with a near-duplicate in arc_length.rs) into solver/sparse_tangent.rs. The shared tangent_free_sparse takes the already-extracted free block (the arc_length signature); corotational's call sites inline the extraction the old private copy did internally. Pure move, no behavior change.
Both contact entry points (2D/3D) dispatch at SPARSE_THRESHOLD to the shared sparse_tangent helpers: CSC conversion, sparse constraint reduction, sparse Cholesky over a fingerprinted symbolic cache. The active set toggles penalty entries per iteration, which changes the CSC pattern — the fingerprinted cache rebuilds the symbolic exactly then and reuses it while the active set is stable. Below the threshold the dense path is byte-identical, error messages included. New parity test: cantilever pressed against a rigid wall through a gap element, coarse mesh (dense) vs fine mesh (sparse) — gap displacement and transmitted force agree within 1e-6.
2D/3D fiber Newton loops dispatch at SPARSE_THRESHOLD to the shared sparse_tangent helpers; modified NR caches the sparse numeric factor of the increment's first tangent (the analogue of the dense cached_l). Dense path below the threshold is unchanged. New parity test: meshed elastic cantilever below/above the threshold, tip deflection vs PL^3/3EI and dense/sparse/modified-NR parity within 1e-4.
The Ernst-iteration solve (2D and 3D) dispatches at SPARSE_THRESHOLD to the shared sparse_tangent helpers with the symbolic Cholesky cached across iterations. New parity test: harp-stayed deck straddling the threshold, cable tension and deck midspan deflection agree within 2%.
The secant-stiffness iteration (2D and 3D) dispatches at SPARSE_THRESHOLD to the shared sparse_tangent helpers. New parity test: soft-clay lateral pile at two mesh densities straddling the threshold (head deflection agrees within 7%; the secant iteration's mesh sensitivity is pre-existing and solver-independent — verified by running the finer mesh through the dense path).
The per-stage solves (2D, 3D, and the staged cable-tension loop) dispatch at SPARSE_THRESHOLD to the shared sparse_tangent helpers. The active-element set changes the sparsity pattern between stages; the fingerprinted symbolic cache rebuilds exactly then. New test: two-stage construction matches a linear reference solve on meshes both sides of the threshold, within 1e-6.
Above SPARSE_THRESHOLD the 2D/3D Newton loops assemble the tangent directly as lower-triangle COO triplets instead of a dense n*n matrix: the element assemblers (truss/frame corotational, springs) are now generic over a scatter closure, so the dense path keeps byte-identical arithmetic. Inclined-support rotation uses the COO triplet helpers from assembly.rs. The final k_dummy dense matrix is only assembled when constraints exist (constraint forces need it); otherwise f_int is accumulated via the triplet drivers. assemble_corotational_public keeps its dense signature (arc-length still consumes it). New test: propped cantilever with an inclined roller straddling the threshold — kinematic restraint check and dense/sparse parity within 1e-4.
The fiber/elastic/spring assemblers are generic over a scatter closure (dense path byte-identical); above SPARSE_THRESHOLD the Newton loops assemble lower-triangle COO triplets and build the CSC free block directly — no dense n*n tangent, no extract_submatrix. tangent_free_sparse_triplets moves to the shared sparse_tangent module (corotational imports it from there now). The constraint-forces tangent rebuild stays dense via the dense closure, mirroring corotational's k_dummy.
… control The sparse path now assembles the co-rotational tangent directly as lower-triangle triplets via assemble_corotational_triplets_2d (now pub(crate)) instead of a dense n*n matrix. Springs go through a scatter closure (dense path byte-identical). The final constraint-forces block uses compute_constraint_forces_sparse over the unreduced CSC free block — no dense K anywhere on the sparse path. Inclined-support behavior is unchanged from the status quo (arc-length never rotated the tangent).
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Sparse-first dynamic analysis paths, on top of the sparse Cholesky series (AMD quotient graph → etree symbolic → supernodal numeric → sparse constraint transform):
a0·M + a1·Kover the union pattern,K_effviaCscMatrix::linear_combination, factored once with sparse Cholesky. Peak reactions keep a sparse full-K (assemble_sparse_2d_ex(build_k_full)) instead of a dense assembly.MAX_DENSE_FALLBACK_DOFS(6000) instead of exhausting the WASM address space.rayleigh_from_modes_sparse).SPARSE_THRESHOLDto shared helpers (solver/sparse_tangent.rs) — CSC conversion, sparse constraint reduction, sparse Cholesky over a fingerprinted symbolic cache reused across Newton iterations (rebuilt only if the pattern changes, e.g. contact active-set or staged element-set changes). Dense paths below the threshold are byte-identical.Tests
main.EqualDOF/RigidLinkconstraints on 2D (sparse branch, nf ≥ SPARSE_THRESHOLD) and 3D; slave==master response parity; dependent-target hard error.linear::solve_2d, Guyan-condensed K == CB K_reduced BB block, slave==master parity.