Skip to content

fix edge case in inline linear solve - #97

Merged
AayushSabharwal merged 2 commits into
mainfrom
os/fix-inline-linear-solve
Jun 10, 2026
Merged

fix edge case in inline linear solve#97
AayushSabharwal merged 2 commits into
mainfrom
os/fix-inline-linear-solve

Conversation

@oscardssmith

Copy link
Copy Markdown
Member

Fix 2 minor bugs found by Fable.

I think this may help with #95

@AayushSabharwal

Copy link
Copy Markdown
Member

Nice finds!

@AayushSabharwal
AayushSabharwal merged commit 87232e8 into main Jun 10, 2026
19 of 23 checks passed
@baggepinnen

Copy link
Copy Markdown
Contributor

Tested this on the model behind #95 (MultibodyComponents.examples.suspension.HalfCar, 22 states, emitted inline linsolve blocks [13, 15, 12, 13, 15, 12, 295]). Test setup: this PR merged on top of my priority-tearing-prototype branch + a probe that marks the joint-axis parameters via maybe_zeros and initializes with LM+pivoted-QR. Three findings:

1. As-is, the PR OOMs on this model

dropzeros!(::SparseMatrixCLIL) calls Base.iszero on the stored values. In get_new_mm the value type here is Num, and iszero(::Num) is a semantic zero test (fraction_iszeroexpand) that explodes on large multibody coefficient expressions. Instrumented run: the process is OOM-killed inside the new dropzeros! at ~row 40 of the 295-row Num-valued mm (DZ! CLIL nrows=295 valtype=Num, then SIGKILL; the pre-PR commit on the same branch completes fine).

Fix that worked for me — a structural-zero hook instead of Base.iszero:

# StateSelection, src/math/sparsematrixclil.jl
cheap_iszero(x) = iszero(x)
# ... in dropzeros!(S::SparseMatrixCLIL):
cheap_iszero(vals[k]) && continue

# ModelingToolkitTearing (has the Symbolics dep)
StateSelection.CLIL.cheap_iszero(x::Num) = SU._iszero(Symbolics.unwrap(x))
StateSelection.CLIL.cheap_iszero(x::SymbolicT) = SU._iszero(x)

Explicit stored zeros produced by sparsevec duplicate-summation are structural Const(0), so the structural check loses nothing.

2. With that fixed, the PR is behavior-neutral on this model

Identical before/after: emitted block sizes ([13,15,12,13,15,12,295]), the 295×295 block (full rank, pivoted-QR relres ≈ 1.4e-14, R-diagonal 1.5e3 → 5e-4), initialization behavior, and the eventual integration failure. So no phantom-zero welds existed here and #95's size problem is unchanged — the fixes look correct (the new unit tests are real bugs), they just don't bite on this model.

3. Incidental finding while measuring (separate issue material)

With a rank-revealing runtime solve in place of \, two of the emitted inline blocks are rank-deficient and inconsistent at runtime, at a fully consistent initialization point (init residual 4e-10):

LSQR n=15 rank=14 relres=0.945   # ‖Ax−b‖/‖b‖ after min-norm pivoted-QR solve
LSQR n=13 rank=12 relres=0.952   # (later, during integration attempts)
LSQR n=13 rank=12 relres=0.0     # rank-deficient but consistent — fine

i.e. the inline-linear-SCC contract "the emitted A\b is solvable at runtime" is violated: b is not in the range of A for the 15×15 block, by ~100%. With the stock \ this produces a SingularException or garbage that silently feeds the RHS; we see immediate Unstable on integration even from a consistent state. Since #94 reworked how eliminated rows contribute to b in __reduce_linear_system!, that's where I'd look first. I'll file this separately with details.

@AayushSabharwal

Copy link
Copy Markdown
Member

Okay wow iszero(::Num) really shouldn't call fraction_iszero

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants