fix: catch the silent symmetric-eigensolver breakdown on a near-singular mass matrix - #140
Conversation
…lar M Both scipy.linalg.eigh and scipy.sparse.linalg.eigsh reduce K x = lam M x through a Cholesky factor of the mass matrix, and that reduction loses accuracy when a very light beam carries a very heavy lump. LAPACK does not raise there, it returns wrong frequencies. A 100 m cantilever with a 4000:1 lump-to-beam ratio reported 0.103 Hz against a true 0.0436 Hz, and the answer wandered non-monotonically with mesh density. solve_modes now measures the backward error of every symmetric solve and redoes it on the general dense path, which factorises neither matrix. Two conditions gate the swap and both are load-bearing. The general result is taken only on a decisive win. A real deck can carry a large backward error without being broken: the bundled NREL 5MW land tower sits at ~2e-2 because its adapter leaves M at cond ~4e10, and there the general path is 1.4x better while splitting a degenerate fore-aft / side-side pair the symmetric solver resolves exactly, which the FA/SS classifier depends on. A genuine breakdown improves by nine orders. Rigid-body modes are excluded from the measurement. For those K x and lam are both ~0, so the relative residual is a ratio of two near-zero quantities and reads ~1 however exact the eigenpair is. Judging a free-free floating solve by the raw maximum condemned it, and acting on that verdict was worse than doing nothing, since the general path filters non-positive eigenvalues and would delete the zero-frequency mode. No existing result changes: the full suite passes with no test edited except one whose workaround for this defect is no longer needed.
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 0d29001144
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| scale = float(np.max(np.abs(eigvals))) | ||
| if scale <= 0.0: | ||
| return 0.0 | ||
| elastic = np.abs(eigvals) > _RIGID_BODY_EIGVAL_RTOL * scale |
There was a problem hiding this comment.
Detect rigid-only subsets without using their own maximum
When n_modes selects only rigid-body modes, their small numerical eigenvalues become the reference scale, so most or all of them satisfy this relative cutoff and are treated as elastic. For example, a symmetric rank-6 stiffness matrix with six free rigid DOFs and n_modes=3 can produce residuals near one, trigger the retry, and switch to _solve_dense_general, which filters non-positive zero modes and returns flexible modes in their place. The rigid-mode test added here avoids this because its requested subset also contains elastic eigenvalues; the cutoff needs an absolute/problem-scale test or must defer the verdict when the returned subset supplies no elastic reference.
Useful? React with 👍 / 👎.
Codex P1 on #140. My rigid-body exclusion keyed its scale off the returned eigenvalues, so a subset containing only rigid-body modes took its own numerical noise as the reference and re-admitted them as elastic, reading their ~1 residuals as a breakdown and retrying into the one path that discards them. Rewriting the classifier to key off strain instead was worse. A genuinely soft mode on a stiff structure carries little strain too: the 0.084 Hz lump mode of a 1e10 N.m^2 beam came in below the cutoff, so the guard stopped firing on the very case it was written for and the answer went back to being 538 percent wrong, silently. Nothing here reliably separates a rigid-body mode from a soft elastic one, so the classification is gone. Both candidate solves are measured the same way over every mode, and the retry still needs a decisive win, so a mode the metric cannot speak to says the same nothing twice and cannot tip the decision. What actually made a false positive dangerous was the fallback dropping zero eigenvalues, changing the mode set rather than just its precision. The retry now runs with keep_rigid_body=True, which clamps rounding-level negatives to zero and keeps them. The asymmetric production path is untouched and keeps the BModes-matching filter it is validated against. A rigid-only subset now returns the modes it was asked for at every size, and the point-mass case is exact again.
|
Fixed, and this one was worth the trouble — the finding was right and chasing it properly turned up a second bug in my own fix. You were correct about the mechanism: a subset containing only rigid-body modes takes its own numerical noise as the reference scale, re-admits them as elastic, and reads their ~1 residuals as a breakdown. I reproduced it — My first fix was your suggested alternative: classify by strain against So I dropped the classification entirely. Nothing here separates a rigid-body mode from a soft elastic one reliably, and both attempts to try made things worse. What made a false positive dangerous in the first place was not the detection but the fallback: it filtered out zero eigenvalues, so a spurious retry changed the mode set rather than just its precision. Fixing that inverts the problem. The retry now runs with The asymmetric production path keeps the BModes-matching Verified: rigid-only subsets return the requested mode count at n = 1, 3, 6 and 10; the point-mass case is exact again (0.00 %); the mesh sweep from 7 to 201 elements converges instead of wandering. 1169 passed, ruff and mypy clean. @codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 7873c56bff
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| gk, gm, n_modes, keep_rigid_body=True, | ||
| ) | ||
| _normalize_columns_l2(alt_vecs) | ||
| alt_worst = _max_residual(gk, gm, alt_vals, alt_vecs) |
There was a problem hiding this comment.
Keep rigid residuals from masking an elastic breakdown
When a near-singular model contains both rigid and elastic modes, taking the maximum across all modes prevents the decisive-win check from detecting corrupted elastic modes: the general candidate still has rigid-mode residuals around 1, so alt_worst cannot be less than 0.1 * worst unless the symmetric residual exceeds roughly 10. Fresh evidence from combining six zero DOFs with the ill-conditioned cantilever and applying an orthogonal basis rotation produced a symmetric fundamental that was 18–44% wrong with elastic residuals of 0.77–0.82, while alt_worst remained about 1 and residual_fallback stayed false. Thus a free-free structure can still silently return the failure this guard is intended to correct; the improvement comparison must avoid allowing unchanged rigid-mode noise to mask a decisive improvement in the elastic modes.
Useful? React with 👍 / 👎.
Codex P1 on #140, and it invalidates the cancellation argument in my last commit. A rigid-body mode reads ~1 in both candidates, so on a maximum it floors the alternative: max(alt_r) stays near 1 and no improvement among the elastic modes can clear a tenth of max(sym_r) unless the symmetric solve is worse than ~10. A free-free model with an elastic mode corrupted to a backward error of ~0.8 sailed straight through, silently, which is the failure the guard exists to catch. The decision is now per mode. Rigid modes contribute ~1 against ~1 and register as no improvement; a corrupted elastic mode contributes ~0.8 against ~1e-9 and registers clearly. Both candidates are sorted over the same spectrum, since the retry already preserves rigid-body modes, so equal indices mean the same mode. An alternative that recovered fewer modes is refused outright: losing one is never an improvement. The regression test builds the case directly, bolting six free rigid DOFs onto the ill-conditioned cantilever and rotating the pair so the blocks are not separable, and pins the mechanism as well as the outcome: the alternative maximum really is floored above 0.1, a maxima comparison really does see no win, and a per-mode one really does.
|
Fixed, and this one directly invalidates the argument I made in the previous round, so thank you for not letting it stand. I claimed the meaningless rigid-body component would "cancel because both candidates are measured the same way". That is true of a per-mode comparison and false of a comparison on maxima, which is what I actually wrote. A rigid mode reads ~1 in both, so it floors The decision is now per mode. Rigid modes contribute ~1 against ~1 and register as no improvement; a corrupted elastic mode contributes ~0.8 against ~1e-9 and registers clearly. Both candidates are sorted over the same spectrum — the retry already preserves rigid-body modes — so equal indices mean the same mode. I also refuse an alternative that recovered fewer modes than the symmetric solve outright, since losing a mode is never an improvement whatever the residuals say. The regression test builds your case directly: six free rigid DOFs bolted onto the ill-conditioned cantilever, the pair rotated by an orthogonal basis so the blocks are not separable by inspection. It pins the mechanism as well as the outcome, asserting that the alternative maximum really is floored above 0.1, that a maxima comparison really does see no decisive win, and that a per-mode one really does — so if anyone later simplifies this back to comparing maxima, that test fails rather than the behaviour silently regressing. 1172 passed, ruff and mypy clean. The three end-to-end behaviours still hold: the point-mass case exact at 0.00 %, rigid-only subsets returning their requested mode count at every size, and the mesh sweep converging. @codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 8eda44ec4a
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| if n_alt < n_sym: | ||
| return np.zeros(0, dtype=bool) |
There was a problem hiding this comment.
Preserve negative modes before comparing by index
When geometric softening makes K indefinite and the caller requests only a subset, _solve_dense_general(..., keep_rigid_body=True) still filters genuinely negative eigenvalues below its rounding floor, then fills the subset with later positive modes. The alternative can therefore have the same length as the symmetric result while describing a shifted spectrum; this check accepts it and the per-index residual test can replace the result, silently dropping an unstable mode and shifting every subsequent mode. Verify eigenvalue correspondence or retain negative modes on the retry rather than treating equal counts as proof that indices match.
Useful? React with 👍 / 👎.
Two things. Codex P2 on #140: preserving only rounding-level zeros still filtered genuinely negative eigenvalues, which an indefinite K produces once run(gravity=...) loads a column past its buckling weight. The alternative then came back the same length, backfilled from higher up, while describing a shifted spectrum, and the per-index comparison was reading two different spectra against each other. eigh filters nothing, so the retry now filters nothing either: keep every real finite eigenvalue, sign included. Both paths return the n_modes smallest real eigenvalues, so equal indices describe the same mode by construction rather than by assumption. The asymmetric production path keeps its positive-only filter. And two of my own tests failed on the Linux BLAS while passing here, because they asserted LAPACK behaviour rather than pyBmodes logic. The degenerate-pair test wanted a split below 1e-9 where one build gives 1e-16 and another 3e-4; it now asks only that the pair is still one mode to engineering precision, which is the property that matters. The masking mechanism is now stated as arithmetic on residual vectors instead of being coaxed out of a matrix pair, so it pins the two decision rules exactly and cannot drift with the vendor library.
|
Fixed, and you were right that equal counts prove nothing. Preserving only rounding-level zeros still filtered genuinely negative eigenvalues, which is not hypothetical here: Rather than verify correspondence after the fact I removed the mismatch at the source. Regression test builds the indefinite case directly and asserts the trap explicitly — the default filter and the preserving one return the same length while the preserving one alone contains the unstable mode, and Separately, CI caught two of my own tests being wrong rather than the code: they asserted LAPACK behaviour and passed on Windows while failing on the Linux BLAS. The degenerate-pair test demanded a split below 1e-9, where one build gives 1e-16 and another 3e-4 — it now asks only that the pair is still one mode to engineering precision, which is the property that actually matters. And the masking mechanism from the last round is now stated as arithmetic on residual vectors rather than coaxed out of a matrix pair, so it pins the two decision rules exactly and cannot drift with the vendor library. 1176 passed locally, ruff and mypy clean. @codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: fbf9e06a84
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| # stays exact there. | ||
| residual_fallback = False | ||
| if sym: | ||
| sym_r = _modal_residuals(gk, gm, eigvals, eigvecs) |
There was a problem hiding this comment.
Measure residuals against the symmetrized matrices
For matrices with asymmetry below symmetry_rtol, the symmetric paths solve 0.5 * (K + K.T) and 0.5 * (M + M.T), but this evaluates their residuals against the original matrices. In a high-dynamic-range model, an accepted skew of 1e-12 * max|K| can be large relative to a soft mode, so an accurate eigh result appears to have residual greater than 0.1; the retry then decisively wins only because eig solves the different, unsymmetrized problem and replaces the spectrum. For example, an accepted 3×3 stiffness with low eigenvalues around 1e-12 and skew 8e-13 shifts the first eigenvalue by about 8–10%. Compute this check against the same symmetrized matrices used by the selected symmetric solver.
Useful? React with 👍 / 👎.
Codex P1 on #140. The symmetric paths symmetrise internally, but the guard judged their modes against the raw matrices. The skew that _is_effectively_symmetric tolerates is only small relative to max|K|, so in a model with a wide dynamic range it can be the same size as a soft mode own eigenvalue. An exact symmetric solve then reads as broken, and eig on those same raw matrices wins decisively purely by answering a different question, replacing the symmetric spectrum the caller was promised with the skewed one. Both the measurement and the retry now use the symmetrised pair. The regression test builds a 3x3 with eigenvalues from 1 down to 1e-12 and skew just inside the tolerance, coupling the two softest modes so the perturbation is comparable to the eigenvalue rather than lost against it. It asserts the pair is accepted as symmetric, that no retry fires, that the returned spectrum is the symmetric one, and that measuring the very same modes against the raw matrices would have read above the threshold.
|
Fixed. This one was the subtlest yet and the reasoning is exactly right. The symmetric paths symmetrise internally, and I was judging their output against the raw matrices — so the guard was asking a question the solver had never been asked. The skew Both the measurement and the retry now use the symmetrised pair, so the whole comparison is about one problem. The regression test builds your construction: a 3×3 with eigenvalues from 1 down to 1e-12 and skew just inside the tolerance. Two details I had to get right and mention in case they matter to the next reader — the symmetry test compares 1179 passed, ruff and mypy clean. @codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: b724726aeb
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| gk_s = 0.5 * (gk + gk.T) | ||
| gm_s = 0.5 * (gm + gm.T) | ||
| sym_r = _modal_residuals(gk_s, gm_s, eigvals, eigvecs) |
There was a problem hiding this comment.
Measure diagnostics against the symmetrized pair
When accepted skew is present, these symmetrized matrices correctly drive the retry decision, but _build_diagnostics is later passed the raw gk/gm and recomputes residuals against a problem the returned modes did not solve. With the added 3×3 accepted-skew construction, the correct symmetric result reports residual_fallback=False while diag.max_residual is about 0.37, falsely flagging a healthy solve as defective in the certification telemetry. Reuse the residuals from the solved pair or pass the symmetrized matrices when building diagnostics.
Useful? React with 👍 / 👎.
Codex P2 on #140, and a direct consequence of a call I made in the last commit. I moved the retry decision onto the symmetrised matrices but deliberately left _build_diagnostics on the raw ones, to avoid perturbing a documented field. That was the wrong trade: with accepted skew present, a correct symmetric solve reported max_residual around 0.37, so the field advertised as certification telemetry was flagging a healthy result as defective. solve_modes now names the pair the returned modes actually solve once, symmetrised on the symmetric paths and raw on the general one, and uses it for the retry decision and the diagnostics alike. The field docstring says which basis it is on, since a backward error is meaningless without naming the problem it is against. Genuinely symmetric decks are unaffected: with no skew the symmetrised pair is the raw pair.
|
Fixed, and this one is squarely on me rather than an oversight — I considered exactly this last round and decided against it. My reasoning was "don't perturb a documented field", so I moved the retry decision onto the symmetrised matrices and deliberately left
Genuinely symmetric decks are unaffected, since with no skew the symmetrised pair is the raw pair, which is why nothing in the existing suite moved. 1180 passed, ruff and mypy clean, strict docs build clean. Six rounds, six real defects, every one in code I had convinced myself was right — including this one, which I had actively reasoned about and got backwards. Thank you for the persistence; this guard is materially better for it. @codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 55e30ef985
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| & np.isfinite(eigvals_real.real) | ||
| & (eigvals_real.real > 0.0) | ||
| ) | ||
| valid = np.isreal(eigvals_real) & np.isfinite(eigvals_real.real) |
There was a problem hiding this comment.
Preserve numerically complex real-spectrum modes on retry
When retrying a highly ill-conditioned symmetric free-free problem, QZ can represent theoretically real zero modes as small complex-conjugate pairs that np.real_if_close(..., tol=1000) does not coerce, so this mask still deletes them even with preserve_full_spectrum=True. Because enough later real eigenvalues remain to fill n_modes, the equal-count guard passes and elastic modes backfill the removed indices before the fallback is accepted. Fresh evidence from the current TestRigidModesCannotMaskAnElasticBreakdown._rigid_plus_ill_conditioned fixture is that requesting 10 modes returns only 2 of its 6 rigid modes on this build, silently shifting the spectrum despite the new preservation path.
Useful? React with 👍 / 👎.
Codex P1 on #140, which I had independently reached from the other direction while reviewing my own work. QZ can represent a symmetric problem theoretically real eigenvalues as small complex-conjugate pairs that real_if_close will not coerce, so the real filter deletes them even with no sign filter, and a truncated request backfills the gap from higher up. Equal counts then prove nothing about equal indices, which was the assumption the whole per-index comparison rested on. My first attempt demanded that every mode survive. That is too strict to be useful: on this machine eig drops six of 249 on a healthy fixture, all at the stiff end, far above anything a caller asks for. The invariant that actually matters is narrower, so the retry now verifies it directly. Nothing discarded may fall inside the returned window, and when something does the retry declines rather than swapping. Codex evidence shows why declining is the right answer rather than a weaker one: on their build the dropped pairs are the zero modes themselves, and keeping the alternative would backfill them with elastic modes and silently shift the spectrum. A guard added to stop a silent wrong answer must not be able to introduce one. That leaves the guard reliable and platform-independent for the case it was built for, a near-singular mass matrix with no rigid-body modes, and safe but not always effective when rigid-body modes coincide with one. Said plainly in the module docstring, the CHANGELOG and the test that now asserts the portable property there rather than a build-dependent rescue. The retry also gets its own entry point rather than a flag on _solve_dense_general, so the asymmetric production path is provably untouched.
|
Confirmed and fixed — and we arrived at this one from opposite directions simultaneously. I was auditing my own retry for exactly this hole when your comment landed; your evidence is the more alarming half. I had reproduced the mechanism but drawn the wrong conclusion from it. On this machine The invariant is narrower than either: nothing discarded may fall inside the returned window. The retry now checks that directly and declines to swap when it cannot verify it. Declining rather than reaching for something cleverer is deliberate, and your evidence is the argument for it — backfilling a dropped zero mode with an elastic one would be precisely the silent spectrum shift this guard exists to prevent, and a guard added to stop a silent wrong answer must not be able to introduce one. That leaves an honest scope statement rather than a claim of generality, now written into the module docstring, the CHANGELOG, and the tests:
The rigid+ill-conditioned test now asserts the portable property — the requested modes come back, ascending, none lost — instead of a rescue that is genuinely build-dependent. The rescue itself stays pinned on the unambiguous case. The retry also gets its own entry point instead of a flag on 1183 passed, ruff and mypy clean. @codex review |
The drop-above-the-window test was reading a real matrix and trusting that its LAPACK build would emit the complex pair at the stiff end. Which end those land on is exactly what differs between builds, and is what CI caught me on twice already. An explicit rotation block places the pair where the test needs it and says the same thing everywhere.
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 0fc9e0cb94
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| _decisively_improved_modes( | ||
| sym_r, alt_r, alt_vals.size, eigvals.size, | ||
| ) |
There was a problem hiding this comment.
Match the sparse spectrum before comparing residuals
When a symmetric system takes the sparse path and has negative modes, eigsh(..., sigma=0, which="LM") returns modes nearest zero, whereas _general_spectrum_for_retry sorts algebraically and returns the most negative modes first. If the sparse residual exceeds the threshold—such as for a large, near-singular, post-buckling model—the per-index comparison therefore pairs different modes, and an exact far-negative alternative can trigger a swap even when the corresponding sparse modes did not improve. Fresh evidence beyond the earlier sign-filter issue is that the current retry also runs after sparse_shift_invert, whose selection window differs from the retry's window; verify eigenvalue correspondence or make both paths select the same spectrum before using this mask.
Useful? React with 👍 / 👎.
| alt_vals, alt_vecs, ordering_sound = _general_spectrum_for_retry( | ||
| gk_s, gm_s, n_modes, | ||
| ) |
There was a problem hiding this comment.
Preserve the symmetric result when the retry fails
When the symmetric solve returns a residual above the threshold but scipy.linalg.eig fails to converge on the same ill-conditioned pencil, this unguarded retry raises LinAlgError and discards the already available symmetric result. That turns the guard from a safe, best-effort correction into a new hard failure precisely on the defective inputs most likely to make the alternative solver fail; catch retry failures and decline the swap while retaining the original result and diagnostics.
Useful? React with 👍 / 👎.
Two from Codex on #140. P1, the sparse spectrum. eigsh(sigma=0, which="LM") selects the modes nearest zero in magnitude while the retry selects the algebraically smallest, so with negative eigenvalues present, a post-buckling gravity column, the two windows are different sets and a per-index comparison pairs unrelated modes. Rather than reconcile the two selections I removed the case. Only the dense symmetric path is retried now, which is a statement about which matrix each routine factorises rather than a convenience: eigh reduces through a Cholesky factor of the mass matrix, the failure this guard exists for, while eigsh factorises K and is unaffected. The mesh sweep that motivated the whole change already showed that, returning correct frequencies on exactly the meshes large enough to take the sparse path. So the mismatch is gone by construction and nothing is lost. P2, the retry could raise. A pencil defective enough to break the symmetric reduction can also break eig, and an unguarded LinAlgError would have made this guard destroy usable results on precisely the inputs it was added to help. It now declines and keeps the symmetric result, logging why.
|
Both fixed. The first one had a root-cause fix that removes the class rather than handling it. Sparse spectrum mismatch (P1). Correct, and the diagnosis is exact: Rather than reconcile the two selections, I removed the case. The retry is now scoped to the dense symmetric path only — and that is a statement about which matrix each routine factorises rather than a convenience. Unguarded retry (P2). Also correct, and the framing is the right one: a guard added to prevent a silent wrong answer must not be able to convert a usable result into a hard failure, least of all on the defective inputs most likely to break the second solver too. The retry now catches 1187 passed, ruff and mypy clean. @codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: fde9191976
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| stacklevel=2, | ||
| ) | ||
| eigvals, eigvecs = alt_vals, alt_vecs | ||
| path = "dense_general" |
There was a problem hiding this comment.
Avoid reporting an overlarge request as filtered modes
When a dense symmetric problem accepts the residual retry and n_modes exceeds the matrix dimension, setting path to dense_general makes the later mode-count warning claim that the shortfall came from complex or non-positive eigenvalues. The retry helper actually preserves non-positive eigenvalues, and the shortfall is simply the benign min(n_modes, ngd) truncation that the surrounding warning logic intends to exclude; for example, a 117-DOF ill-conditioned system requested with 1000 modes returns all 117 but reports the pencil as defective. Gate that warning on modes actually discarded by the general solver rather than the post-retry path label.
Useful? React with 👍 / 👎.
| `solve_modes` now checks the backward error `||K x - λ M x|| / ||K x||` | ||
| of every symmetric solve and, when it is large, redoes it through the | ||
| general dense path, which factorises neither matrix. The retried result |
There was a problem hiding this comment.
Document the dense-only retry scope
The release note says every symmetric solve is checked and redone through the general path, but the implementation now deliberately restricts retries to path == "dense_symmetric"; sparse results are never retried even when their residual exceeds the threshold. This also conflicts with the preceding claim that eigsh factorizes the mass matrix, whereas the implemented rationale says shift-invert factorizes K. Update the public change description so users do not expect residual_fallback or a corrective retry on sparse solves.
Useful? React with 👍 / 👎.
A self-review pass plus two from Codex on #140. Found while auditing: the retry was reachable at any size. A sparse solve that fails to converge falls back to dense however large the system, and an unbounded dense eig there could spend minutes on a result already in hand. A guard against a silent wrong answer should not be able to turn one into a silent hang, so it now declines above residual_retry_max_ndof and leaves the backward error to the diagnostics. Codex, mode-count warning. After a swap the path label reads dense_general, which made the shortfall warning claim modes had been filtered as complex or non-positive. The retry preserves the whole spectrum, so nothing was filtered, and a 117-DOF system asked for 1000 modes was reported as defective for returning its 117. The warning now keys on modes actually discarded rather than on the label, which also fixes the same spurious report on the pre-existing general path. Codex, CHANGELOG. The note still claimed every symmetric solve is retried and that eigsh factorises the mass matrix, both of which my own scoping change had made false. Rewritten to say the retry is dense-only, that sparse results never set residual_fallback, and why. Also states in the diagnostics that symmetric=True alongside path=dense_general is not a contradiction: the first describes the input, the second which routine ran, and residual_fallback separates that case from a genuinely asymmetric solve.
|
Both fixed, plus one I found in a self-review pass before it was reported. Mode-count warning (P2). Correct. After a swap the path label reads CHANGELOG (P2). Also correct, and a good catch on a doc claim I had invalidated myself: my own scoping change made "every symmetric solve" and "eigsh factorises the mass matrix" both false, and I updated the code without updating the note. Rewritten to say the retry is dense-only, that sparse results never set Found in self-review: the retry was unbounded in size. A sparse solve that fails to converge falls back to dense at any size, so the retry was reachable on a 5000-DOF system where a dense I also ran an adversarial sweep over the guard's edge cases: 1192 passed, ruff and mypy clean. @codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 603eae8b50
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| if ordering_sound | ||
| else np.zeros(0, dtype=bool) | ||
| ) | ||
| if improved.any(): |
There was a problem hiding this comment.
Reject retries that degrade other requested modes
When QZ decisively improves one mode but makes another requested mode less accurate, improved.any() accepts the retry and lines 419–421 replace the entire spectrum. For example, on a 4×4 SPD pencil with cond(M) ≈ 1e16, the symmetric residuals can be [0.52, 1.1e-3, 5.7e-6, 0.017] while the general residuals are [1.1e-10, 3.0e-11, 3.5e-5, 0.26]; this condition swaps because mode 0 improved, thereby replacing a previously acceptable fourth mode with one above the failure threshold. Require the whole candidate spectrum to be non-regressive, rather than accepting it when any single mode improves.
Useful? React with 👍 / 👎.
Codex P1. Acceptance keyed on "any mode improved decisively", but taking the retry replaces the whole spectrum, not the modes that prompted it. So a candidate that rescued mode 0 while pushing a previously acceptable mode 3 above the failure threshold was accepted, handing back a new bad mode in place of an old one. The comparison now returns both verdicts and the caller requires improvement somewhere and regression nowhere. Regression deliberately mirrors improvement rather than introducing a second notion of acceptable: a mode has regressed when it ends up above the threshold and is worse there by the margin that would have counted as decisive the other way. The symmetry is what keeps rigid-body modes out of it, since their residual reads ~1 in both candidates and wobbles either way, and a bare alt > sym test would read that as a regression and block every rescue sitting beside a free-free mode. Tests take the exact residual vectors from the report, plus the two cases the rule has to keep apart: a mode six times worse but nowhere near the threshold is not a regression, and neither is rigid-body noise.
|
Fixed, and this is the sharpest finding of the series — it catches a gap in the shape of the rule rather than in its plumbing. The rule keyed on "any mode improved decisively", but taking the retry replaces the whole spectrum, not the modes that prompted it. So a candidate that rescued mode 0 while pushing a previously acceptable mode 3 above the failure threshold was accepted, handing back a new bad mode in place of an old one. Your residual vectors are now the test. The comparison returns both verdicts and the caller requires improvement somewhere and regression nowhere. One design note on how I defined regression, since the obvious version is wrong. A bare
Both are pinned, along with the rigid-body noise case and an end-to-end test that corrupts one mode of the candidate and asserts the swap is declined. 1196 passed, ruff and mypy clean. @codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: dc26e04879
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| factor = _SOLVER_OPTIONS.residual_retry_improvement | ||
| sym, alt = sym_r[:n], alt_r[:n] | ||
| improved = (sym > threshold) & (alt < factor * sym) | ||
| regressed = (alt > threshold) & (sym < factor * alt) |
There was a problem hiding this comment.
Reject candidates that cross the residual threshold
When a previously acceptable mode crosses the threshold by less than the configured tenfold factor, this condition does not mark it as regressed, so an improvement elsewhere still replaces the whole spectrum with a newly failing mode. For example, sym_r=[0.52, 0.09] and alt_r=[1e-10, 0.8] produce an improvement and no regression, even though the second mode moves from below the 0.1 threshold to a residual of 0.8; this also does not involve the rigid-body-noise case because those modes already start above the threshold. Treat an acceptable-to-failing threshold crossing as a regression in addition to the decisive-worsening comparison.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 430ec073dc
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| # retry decision below — have to be measured against that pair, not | ||
| # against the raw one. Reporting the raw backward error would flag a | ||
| # correct solve as defective in telemetry meant to be auditable. | ||
| res_k, res_m = (0.5 * (gk + gk.T), 0.5 * (gm + gm.T)) if sym else (gk, gm) |
There was a problem hiding this comment.
Avoid rematerializing matrices after ordinary sparse solves
When a large symmetric problem successfully takes sparse_shift_invert with return_diagnostics=False, this unconditionally constructs two full dense symmetrized matrices even though the retry is excluded by path == "dense_symmetric" and diagnostics are never built. _solve_sparse_shift_invert has already symmetrized the same pair internally, so every ordinary sparse solve now pays for a second pair of O(ngd²) allocations and matrix additions, adding substantial memory traffic on exactly the large systems for which the sparse path exists. Construct res_k/res_m lazily only when the dense retry or returned diagnostics needs them.
Useful? React with 👍 / 👎.
Codex P2. Measuring against the symmetrised pair was done by building it, unconditionally, on every solve. A large sparse solve therefore paid for two dense ngd-square allocations it never used, on exactly the systems the sparse path exists to keep cheap. I noted this cost in my own audit and judged it acceptable, which was wrong: at ngd = 1500 it is a 54 MB peak against 0.3 MB, a factor of 187. 0.5 (A + A.T) v equals 0.5 (A v + A.T v), so the basis change is now carried as a flag and applied through two thin products. Identical results to 1e-12, pinned by a test, along with one asserting the flag is not inert and one bounding the allocation below a quarter of a single matrix. The symmetrised pair is still built in one place, the retry, because eig needs matrices rather than products. That branch is rare and already size-capped.
|
Fixed, and this one I had already found and misjudged — which makes it the more useful of the two outcomes. During my hardening pass I noted exactly this ("two O(n²) matrix adds on every solve… acceptable but wasteful when the retry never runs") and moved on. It is not acceptable: at ngd = 1500 the eager symmetrisation peaks at 54 MB against 0.3 MB, a factor of 187, and it falls on precisely the large sparse solves the sparse path exists to keep cheap. Worse, I introduced it — before this PR the diagnostics measured against the raw matrices and allocated nothing. The identity The symmetrised pair is still materialised in one place — the retry itself, because Worth noting what this says about my hardening pass. I ran it specifically to pre-empt findings, and it did surface this one — I just weighed it wrong, because I was thinking about the retry branch where the cost is bounded and rare, rather than the diagnostics path where it is unconditional. Noticing a cost and mis-scoping who pays it is a distinct failure from not noticing it, and measuring rather than estimating would have settled it in a minute. (Reposted: my previous comment lost two code spans to shell substitution.) 1219 passed, ruff and mypy clean. @codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: acb57844de
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| gk_s = 0.5 * (gk + gk.T) | ||
| gm_s = 0.5 * (gm + gm.T) |
There was a problem hiding this comment.
Delay symmetrization until the retry actually runs
For every eligible dense symmetric solve—including healthy solves and large sparse-to-dense fallbacks up to 2000 DOFs—these lines still allocate both full symmetrized matrices before the residual is compared with the retry threshold. Fresh evidence beyond the earlier sparse-path finding is that the new thin-product helper is not used for this initial check, so a retry that never runs still pays the repeated O(ngd²) memory traffic; compute sym_r from the raw pair with symmetrise=True, then materialize gk_s and gm_s only inside the threshold-triggered branch where eig needs them.
Useful? React with 👍 / 👎.
Codex P2, and the same mistake as the previous commit one level in. I switched the diagnostics to thin products and left the sibling call site a hundred lines above still building the pair, before the threshold that decides whether the retry runs at all has even been tested. Every healthy dense solve up to the size cap paid two dense allocations to find out it did not need them. The pre-threshold check now goes through the products too. The pair is built inside the branch, where eig genuinely needs matrices. The test I first wrote for this asserted the whole solve allocates less than half a matrix, which is impossible and would have been a bad test even if it passed: eigh takes matrices, so _solve_dense_symmetric legitimately builds the pair it needs, and that dominates any peak measurement. It now asserts by identity that the check reads the caller arrays rather than a copy, which is the property actually at stake.
|
Fixed — and this is the same mistake as the previous commit, one level in. I switched the diagnostics to thin products and left the sibling call site a hundred lines above still building the pair, before the threshold that decides whether the retry runs at all has even been tested. So every healthy dense solve up to the size cap paid two dense allocations to discover it did not need them. The pre-threshold check now goes through the products too; the pair is built inside the branch, where Two things worth recording, since both are about how I verified rather than what I changed. The test I first wrote for this was wrong, and would have been bad even if it passed. I asserted the whole solve allocates less than half a matrix. That is impossible — My previous fix was partial in a way I should have caught. I applied the helper at one call site and did not check the others. Grepping for the pattern I had just introduced would have found this in seconds, and I did not do it. 1220 passed, ruff and mypy clean. @codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: dfbb96c99d
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| # Through the products, not the pair: this runs on every eligible | ||
| # solve, healthy ones included, and the threshold below has not | ||
| # been tested yet. | ||
| sym_r = _modal_residuals(gk, gm, eigvals, eigvecs, symmetrise=True) |
There was a problem hiding this comment.
Avoid full-width products for full-spectrum solves
When n_modes is omitted—the public default—or requests most of the spectrum, eigvecs is ngd × ngd, so these are not thin products: _apply evaluates both A @ V and A.T @ V, creating several full matrix-sized temporaries for each of K and M. Thus every healthy dense solve within the 2000-DOF retry cap still incurs substantial peak memory and four extra dense matrix multiplications before the threshold is checked; locally, the full-width residual call peaked at about five times one matrix allocation, whereas the added test covers only five vectors. Handle the full-width case separately or reuse the symmetrized matrices/products already required by the dense solve.
Useful? React with 👍 / 👎.
The two forms of sym(A) @ V peak at 3nk and 2n^2 bytes of temporaries. Splitting the product was only ever cheaper for a narrow block, but n_modes=None is the public default and returns the whole spectrum, so the wide case was paying for two full matmuls to avoid one copy. Select the route by width. The crossover is measured rather than derived: a flop-count estimate put it at k = n/3 when the allocator puts it at 2n/3, which would have taken the dearer route across a third of the range.
|
Fixed in e93347f — and the measurement corrected my first attempt, so it's worth recording what it said. You're right that the products stop being thin. The two routes peak at What I got wrong on the first pass was where to put the switch. A flop-count estimate said The real crossover is Two tests cover it. One asserts I did consider reusing the symmetrised matrices from the dense solve. The residual measurement runs before we know whether a retry is warranted, and on the healthy path there is no retry, so materialising there would put the cost back on every solve to serve the rare one. Choosing by width keeps the common case cheap without a second code path. |
The width test compared three live tracemalloc peaks. That measures what the allocator did, not what the routes cost: a block freed by an earlier call gets reused by a later one, so the three time-share buffers differently depending on measurement order. It passed on Windows and failed on Linux at k/n = 0.9, where the route taken was correct and only the comparison figure had shifted. Lift the decision into _prefer_materialised and assert it against the cost model, which is deterministic everywhere. Adds explicit cover for the full-width case that prompted the rule and for the thin case it must not disturb, and widens the numerical-agreement test to straddle the crossover.
|
CI caught the width test on Linux while it passed on Windows, so that test is gone in cad2de1 — and the failure is worth recording, because the route it flagged was correct. It compared three live Same trap as the two LAPACK-dependent tests earlier in this PR: asserting measured platform behaviour rather than the thing under test. So the decision is now a pure predicate, The measurement that established the model stays in the comment above the predicate and in my previous reply, which is the right place for it — it justified the constant, it should not be re-run as an assertion. |
The refusal is on a mode getting strictly worse, so an accepted retry proves an acceptable mode did not get worse — not that it improved. Below the regression floor it may move either way, which is what test_a_mode_that_worsens_but_stays_acceptable_is_not_a_regression pins. Say that rather than the slightly stronger claim.
|
@codex review |
|
@codex review Four commits since the last review, all responding to the full-width residual finding:
CI is green on |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 96040ec551
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| if v.ndim > 1 and _prefer_materialised(a.shape[0], v.shape[1]): | ||
| return np.asarray(0.5 * (a + a.T) @ v) | ||
| return np.asarray(0.5 * (a @ v + a.T @ v)) |
There was a problem hiding this comment.
Reduce the full-width residual peak, not only the product peak
When n_modes=None or most modes are requested, this routing reduces the number of full matrix multiplications but does not resolve the reported peak-memory problem. Fresh evidence in the new route is that its 3nk versus 2n² model stops at _apply: _modal_residuals retains the full-width kx while constructing mx, then retains both while allocating mx * eigvals and the subtraction result, so the call still simultaneously holds roughly four additional n×n arrays (over 100 MB at the 2000-DOF retry cap, excluding the eigenvectors). Compute these residuals in blocks or reuse/in-place the products so the public full-spectrum default no longer incurs that peak.
Useful? React with 👍 / 👎.
Routing the product form by width bounded _apply but not its caller. _modal_residuals held kx while building mx, then both while forming mx * eigvals and the difference, so a full-spectrum request kept four ngd-square arrays live purely to measure: 128 MB at the 2000-DOF retry cap, measured. The residual is per mode, so no step needs every mode at once. Sweep in column blocks and accumulate the two norm vectors. Peak is then set by the block width rather than by an argument the caller chooses, and the thin case stays a single pass computed exactly as before. At ngd = 2000 with every mode requested the peak drops from 128 MB to 10.3 MB, results agreeing to 1e-12. The same measurement shows no change at all for six modes. This also removes the need for the width test in _apply, since blocking guarantees it never sees a wide block.
|
Confirmed and fixed in 022b06c. The finding was exactly right, including the figure: I measured 128 MB at the 2000-DOF cap. My previous fix bounded The way out is that the residual is per mode, so no step needs every mode present at once. The sweep now runs in column blocks and accumulates the two norm vectors, which ties the peak to the block width instead of to an argument the caller chooses. Two things worth drawing out. The common thin case is completely untouched — same peak, same time, one pass per matrix, computed exactly as before. And the full-width case costs 0.13 s more, because blocked products do This also makes the width test I added last round unnecessary, so Tests follow the same principle as last round: the invariant is asserted by spying on the width every pass actually sees, not by comparing I also A/B'd the solver against |
128 was a round number. Measured across widths on a 2000-DOF full spectrum it turns out to sit at the knee: 1.4 s at 16 and 32 columns where per-call BLAS overhead dominates, 0.76 s at 64, 0.42 s at 128, against a 0.38 s floor that 256 and above buy with two to twelve times the peak. Record that where the constant is defined. Add a test that the width changes nothing but the peak, so it can be retuned freely. The residual is per mode, so blocking only partitions independent columns.
|
@codex review |
|
@codex review Two commits since the last pass, both on the residual sweep.
CI is green on |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 4dd48df99f
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| # 2000-DOF system took 1.4 s at 16 and 32 columns, where per-call BLAS | ||
| # overhead dominates, then 0.76 s at 64 and 0.42 s at 128, against a | ||
| # 0.38 s floor that 256 and above buy with two to twelve times the peak. | ||
| _RESIDUAL_BLOCK = 128 |
There was a problem hiding this comment.
Scale the residual block to small matrices
When ngd < 192 and the caller uses the public n_modes=None default—or otherwise requests most modes—the first block can exceed the 2n/3 crossover documented immediately above. Thus removing _prefer_materialised does not actually guarantee that _apply receives a narrow block: for a 100-DOF full-spectrum solve it receives a square block and performs two GEMMs per matrix instead of the cheaper materialized-symmetry route, regressing the common small dense case. Cap the block width relative to ngd or retain width-based routing.
Useful? React with 👍 / 👎.
I claimed blocking made the width test in _apply unnecessary. That is wrong below ngd = 192: the block is capped at 128 columns, which is narrow against a large system but not a small one, so a full-spectrum request there still hands _apply a block past the 2n/3 crossover. The two mechanisms answer different questions. Blocking bounds the peak as the mode count grows; the width test picks the cheaper route within whatever block it gets. Keep both. Capping the block relative to ngd was the other way to close this and is worse at both ends, measured: at ngd = 100 it splits one BLAS call into two, 206 us against 116 us, and at ngd = 1000 it lifts the peak from 5.1 MB to 21.3 MB. Restoring the width test costs nothing at large ngd, where the block is already narrow and it does not fire, and at ngd = 100 runs the sweep in 116 us against 178 us for the same peak.
|
Correct, and the over-claim was mine. Fixed in 31b4d64 by restoring the width test. I said blocking made You offered two remedies and they are not equivalent, so I measured both before choosing. Capping the block relative to Tests cover the specific case you found — a sub-192 Full suite 1251 passed, CI green, and I have re-run the integration suite locally against the upstream decks since the solver change: 127 passed, the six BModes cert cases included. |
|
@codex review |
|
@codex review
Of the two remedies you offered I took the second, having measured both. Capping the block relative to Tests cover the case you found directly, plus one pinning that the two mechanisms are not interchangeable so the claim cannot be made again silently. CI green on |
|
Codex Review: Didn't find any major issues. Nice work! Reviewed commit: ℹ️ About Codex in GitHubYour team has set up Codex to review pull requests in this repo. Reviews are triggered when you
If Codex has suggestions, it will comment; otherwise it will react with 👍. Codex can also answer questions or update the PR. Try commenting "@codex address that feedback". |
Fixes a silent wrong answer in the dense symmetric eigensolver, found while working on #35.
The defect
scipy.linalg.eighreducesK x = λ M xthrough a Cholesky factor of the mass matrix. When a very light beam carries a very heavy lump that matrix goes nearly singular, the reduction loses accuracy, and LAPACK does not raise — it returns wrong frequencies. On a 100 m cantilever at a 4000:1 lump-to-beam ratio the reported fundamental was 0.103 Hz against a true 0.0436 Hz, and the answer wandered non-monotonically with mesh density (0.086 → 0.084 → 0.069 → 0.113 → 0.0857). It predates this branch: it reproduces onmasterwith plaintip_mass.solve_modesnow measures the backward error of a dense symmetric solve and, when it is large, solves again through the general path and compares. It warns when it swaps, andSolverDiagnosticsgainsresidual_fallback.Why the acceptance rule looks the way it does
Every clause is there because a simpler version was demonstrably wrong. This PR went through fifteen review rounds; the rule below is what survived.
eigsh(sigma=0)factorisesK, so it does not have this failure; and itswhich="LM"window is a different mode set that must never be index-comparedeig"win" by answering a different questioneigraises, or the system is too large, or the ordering cannot be verified, the symmetric result standsThree attempts to identify rigid-body modes and exclude them were abandoned: by eigenvalue scale (a rigid-only subset becomes its own reference), by strain (a genuinely soft mode has little of it too — this one silently disabled the guard and returned the 538 %-wrong answer again), and by which side of the failure threshold the value falls on.
What separates the two populations is measured rather than argued: genuine rescues improve by 1e5 to 1e10, rigid roundoff by 11× to 16×. Those numbers are a test, so the constants cannot drift from their evidence.
No existing result changes
The decisive-improvement condition is what guarantees it. On the bundled NREL 5MW land deck — whose adapter leaves
Mat cond ~4e10 — the general path is only 1.4× better while splitting a degenerate fore-aft / side-side pair the symmetric solver resolves exactly, which the FA/SS classifier depends on. The full suite passes with no test edited except one whose workaround for this very defect is no longer needed.What this does not promise
Reliable and platform-independent for the case it was built for: a near-singular mass matrix with no rigid-body modes. Elsewhere safe but not always effective — where rigid modes and a near-singular mass coincide, QZ may return the theoretically real zero modes as complex pairs whose position in the spectrum differs between LAPACK builds, and the retry then declines rather than guessing. Stated in the module docstring, the CHANGELOG and the tests rather than left to be discovered.
Verification
VALIDATION.mdrows: closed-form recovery at a 4e5:1 mass ratio, and the rigid-body non-regression.tests/fem/test_ill_conditioned_mass.pypins each rule and the case that forced it, including that the unguarded symmetric result is more than 25 % out — so the guard cannot be quietly removed.Also fixes a
docs/RELEASE_CHECKLIST.mdreference that moved to the Sphinx tree, and two tests that asserted LAPACK behaviour rather than pyBmodes logic and failed on the Linux BLAS.