diff --git a/CHANGELOG.md b/CHANGELOG.md index 4687948..f55ab4a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,7 +8,100 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). ## [Unreleased] -(nothing yet) +### Fixed + +- **The dense symmetric eigensolver could return confidently wrong low + modes on a near-singular mass matrix, silently.** `scipy.linalg.eigh` + reduces `K x = λ 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. + On a 100 m cantilever with a 4000:1 lump-to-beam mass ratio the + reported fundamental was 0.103 Hz against a true 0.0436 Hz, a factor of + 2.4, and the answer wandered non-monotonically with mesh density. + + `solve_modes` now checks the backward error `||K x - λ M x|| / ||K x||` + of a **dense** symmetric solve and, when it is large, redoes it through + the general dense path, which factorises neither matrix. The retried + result is taken only when it **resolves** a mode the symmetric solve + had failed — improving it by more than 1000× and reaching a backward + error of 1e-3 or better. Both conditions are needed because a + rigid-body mode's residual divides one roundoff quantity by another, + so its value is arbitrary (0.076, 0.79 and 12.4 have all been measured + on healthy models) while its improvement ratio stays near 10x, four + orders short of the 1e5 to 1e10 a real rescue achieves. A `RuntimeWarning` names + the swap, and attributes it to the mass matrix only when the mass + conditioning supports that — a wide stiffness range trips the same + guard with a perfectly conditioned mass. `SolverDiagnostics` gains + `residual_fallback` recording it. + + The **sparse** path is deliberately not retried and never sets + `residual_fallback`. `eigsh(sigma=0, mode='normal')` factorises `K` + rather than the mass matrix, so a near-singular `M` does not degrade + it — on the mesh sweep that motivated this work it returned correct + frequencies on exactly the meshes large enough to select it. Retrying + it would also mean comparing two different mode sets, since + `which="LM"` selects the modes nearest zero in magnitude while the + retry selects the algebraically smallest. + + **No existing result changes.** The decisive-improvement condition is + what guarantees that: a real deck can carry a large backward error + without being broken, and on the bundled NREL 5MW land tower (whose + adapter leaves the mass matrix at cond ~4e10) the general path is only + 1.4× better while *splitting* a degenerate fore-aft / side-side pair + the symmetric solver resolves exactly. + + Acceptance requires the candidate to be **non-regressive** as well as + decisively better somewhere. Taking the retry replaces the whole + spectrum, not the modes that prompted it, so a candidate that rescues + one mode while pushing a previously acceptable one above the failure + threshold is refused — it would hand back a new bad mode in place of an + old one. The guarantee is one-sided and precise: a mode that was + acceptable can end up above the regression floor only by not having got + worse, never as collateral of another mode's rescue. Below that floor + it is free to move either way, which is deliberate — a residual already + that small is not a claim about accuracy worth defending. A mode + already failing carries no verdict either way — above the threshold + neither candidate is trustworthy, and a rigid-body mode, whose residual + divides one roundoff quantity by another and has been measured at 12.4 + against 0.79 on a healthy model, lives entirely in that region. + `max_residual` still reports it. + + The comparison is made **per mode** rather than on the two maxima, so + that rigid-body modes cannot distort it. Their backward error is a + ratio of two near-zero quantities and reads ~1 in both candidates + however exact each is; on a maximum that floors the alternative and + hides a genuinely corrupted elastic mode sitting alongside them, while + per mode it simply registers as no improvement. The retry preserves + zero and negative eigenvalues and verifies that nothing was dropped + from inside the returned window, so it can never backfill a missing + mode with a higher one and shift the spectrum. + + **Scope.** The rescue is reliable and platform-independent for the case + it was built for, a near-singular mass matrix with no rigid-body modes. + Where rigid-body modes and a near-singular mass matrix coincide it is + safe but not always effective: QZ may return the theoretically real + zero modes as complex-conjugate pairs, and where those land differs + between LAPACK builds. When they fall inside the requested window the + alternative's ordering cannot be verified and the retry declines, + leaving the result no worse than before with `max_residual` still + reporting the problem. Declining is deliberate — a guard added to stop + a silent wrong answer must not be able to introduce one. + + The retry is also bounded in size, since a sparse solve that fails to + converge falls back to dense at any size and an unbounded `eig` there + could take minutes on a result already in hand. And it can decline: if + the alternative solver raises on the same defective pencil, the + symmetric result and its diagnostics are kept rather than the whole + solve failing. + +- `SolverOptions` gains five fields for the conditions above, each + governing one of them: `residual_retry_threshold` (what counts as a + failing mode), `residual_retry_improvement` (how much better the + candidate must be for the win to be a rescue rather than roundoff), + `residual_retry_resolved` (the backward error it must actually reach), + `residual_regression_floor` (where a worsened mode has to land before + the worsening counts) and `residual_retry_max_ndof` (the size above + which the retry is not attempted at all). ## [1.18.0] — 2026-08-12 diff --git a/VALIDATION.md b/VALIDATION.md index 2e45755..ff2a0c9 100644 --- a/VALIDATION.md +++ b/VALIDATION.md @@ -102,6 +102,8 @@ metrics: | Discrete mid-span point mass (issue #35) | cantilever carrying one lump at station $a$: $f = \sqrt{3EI/(m a^3)}/2\pi$ (Blevins 1979, Table 8-1) | 1st frequency; mesh-position independence | < 0.5 %; coarse-vs-fine < 0.2 % | (within tol) | [`tests/fem/test_gravity_and_point_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_gravity_and_point_mass.py) | no | | Distributed Winkler soil bed vs the lumped mudline condensation (issue #118) | Psaroudakis et al. (2021) / Yu & Amdahl (2023) Eq. 25 is the exact static condensation of a constant-EI pile on a bed of rate $k = D_P E_{SO}$ | coupled 1st frequency, distributed bed vs lumped springs | < 1 % | 0.4 % (the embedded pile inertia the condensed form drops) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | | Distributed Winkler bed converges on the rigid mudline clamp (issue #118) | limit $E_{SO} \to \infty$; residual compliance scales as the elastic length $(4EI/k)^{1/4}$ | 1st frequency at $E_{SO} \times 10^8$ vs the clamped model | < 1 % | (within tol) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | +| Near-singular mass matrix does not silently corrupt the low modes | cantilever carrying one lump: $f = \sqrt{3EI/(m a^3)}/2\pi$ (Blevins 1979, Table 8-1), at a 4e5:1 lump-to-beam mass ratio | 1st frequency, dense and sparse dispatch sizes | < 0.5 % | (within tol; the unguarded symmetric solve is 137 % out) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no | +| Rigid-body modes are not lost to the solver guard | construction (rank-deficient $K$, well-conditioned $M$; requested subsets of 1, 3, 6 and 10 modes) | zero modes retained, requested mode count returned | exact | (within tol) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no | | Deck-reader `n_nodes` refinement (issue #58) | Euler-Bernoulli closed form, modes 1-3, plus self-convergence $n{=}100$ vs $200$ | bending frequencies | < 1 %; self-convergence < 0.2 % | (within tol) | [`tests/test_refine_mesh.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_refine_mesh.py) | no | **Citations** (full author / year forms used in the table above). diff --git a/src/pybmodes/fem/solver.py b/src/pybmodes/fem/solver.py index 3d04605..8e6a582 100644 --- a/src/pybmodes/fem/solver.py +++ b/src/pybmodes/fem/solver.py @@ -32,7 +32,88 @@ when the sparse path fails to converge (logged as a warning). 3. **Dense general** — ``scipy.linalg.eig`` for genuinely asymmetric systems (offshore decks where the rigid-arm transformation makes - the platform-support block non-symmetric). Matches BModes JJ. + the platform-support block non-symmetric). Matches BModes JJ. Also + the retry path when a symmetric solve comes back with a large + backward error — see below. + +The residual retry +------------------ + +``scipy.linalg.eigh`` reduces ``K x = λ M x`` through a Cholesky factor +of the **mass** matrix, and that reduction degrades once ``M`` is nearly +singular — a very light beam carrying a very heavy lump. The failure is +silent: LAPACK returns confidently wrong low modes rather than raising. +On the case in ``tests/fem/test_ill_conditioned_mass.py`` it reported +0.103 Hz against a true 0.0436 Hz. + +:func:`solve_modes` therefore measures the backward error of a dense +symmetric solve and, above +:attr:`~pybmodes.options.SolverOptions.residual_retry_threshold`, solves +again through the general path and compares. Every other rule below +exists because some simpler version of that comparison was wrong. + +**Only the dense path.** ``eigsh(sigma=0, mode='normal')`` factorises +``K``, not ``M``, so the sparse path does not have this failure and is +never retried. That also avoids comparing two different mode sets: its +``which="LM"`` window selects the modes nearest zero in magnitude while +the retry selects the algebraically smallest. + +**Per mode, not on the maxima.** A rigid-body mode's residual divides one +roundoff quantity by another. Taking maxima lets that noise floor the +candidate's worst value and hide a genuinely corrupted elastic mode +beside it. + +**Judged by the size of the win.** No absolute bar separates a rescue +from rigid noise, because the noise value is arbitrary — 0.076, 0.79 and +12.4 have all been measured on healthy models, and the first is *below* +the failure threshold. Identifying such modes was tried three times and +abandoned: not by eigenvalue scale, which a rigid-only subset makes its +own reference; not by strain, which a genuinely soft mode also has +little of; and not by which side of the threshold the value falls on. +What does separate them is the ratio. Rescues improve by 1e5 to 1e10, +roundoff by 11x to 16x, so acceptance needs +:attr:`~pybmodes.options.SolverOptions.residual_retry_improvement` *and* +a candidate that reaches +:attr:`~pybmodes.options.SolverOptions.residual_retry_resolved`. + +**Non-regressive.** Accepting replaces the whole spectrum, so a candidate +that rescues one mode while pushing another past the threshold is a +trade, not an improvement. + +**The same matrices throughout.** Both symmetric paths symmetrise +internally, so the measurement and the retry use the symmetrised pair. +The tolerated skew is small only relative to ``max|K|``, which in a +wide-dynamic-range model can still swamp a soft mode's own eigenvalue; +judging an exact solve against the raw matrices reads as a failure, and +``eig`` on those same matrices then "wins" by answering a different +question. + +**The same spectrum throughout.** The retry keeps every real eigenvalue, +zeros and negatives included, and verifies that nothing was discarded +from inside the returned window. ``eigh`` filters nothing, so any filter +here would return a different set of the same length, backfilled from +higher up, and equal indices would stop meaning equal modes. Both +omissions are reachable: a free-free model's zero modes, and the negative +eigenvalues an indefinite ``K`` produces once ``run(gravity=...)`` loads +a column past its buckling weight. + +**It can always decline.** If the alternative raises on the same +defective pencil, or the system is larger than +:attr:`~pybmodes.options.SolverOptions.residual_retry_max_ndof`, or its +ordering cannot be verified, the symmetric result stands. + +What this does not promise +-------------------------- + +The rescue is reliable and platform-independent for the case it was built +for: a near-singular mass matrix with no rigid-body modes. Elsewhere it +is *safe* but not always *effective*. Where rigid-body modes and a +near-singular mass coincide, QZ may return the theoretically real zero +modes as complex-conjugate pairs, and where those land differs between +LAPACK builds; inside the requested window the ordering cannot be +verified and the retry declines. Declining is deliberate — a guard added +to stop a silent wrong answer must not be able to introduce one — and +``max_residual`` still reports the problem. Note on the user-spec mode choice: ``eigsh(..., sigma=0, mode='buckling')`` reduces to ``OP = K^-1 K = I`` for ``sigma=0``, @@ -82,9 +163,16 @@ class SolverDiagnostics: path : which solver path produced the result. One of ``"sparse_shift_invert"``, ``"dense_symmetric"``, ``"dense_general"``. - symmetric : whether the assembled matrices were treated as symmetric - (``eigh`` / sparse) rather than routed through the general - ``eig`` path. + symmetric : whether the assembled matrices were **classified** as + symmetric, i.e. whether their asymmetry was within + :attr:`~pybmodes.options.SolverOptions.symmetry_rtol`. This is a + property of the input, not a record of which routine ran, so a + residual retry leaves it ``True`` while moving ``path`` to + ``"dense_general"``. That pairing is not a contradiction: the + matrices were symmetric, and the general routine was used on their + symmetrised form because the symmetric one had failed on it. + ``residual_fallback`` is what distinguishes that case from a + genuinely asymmetric solve. n_requested : modes asked for (``None`` means the full spectrum). n_returned : modes actually returned. Fewer than ``n_requested`` means the general path filtered out complex / non-positive @@ -94,12 +182,29 @@ class SolverDiagnostics: attempted and failed, so the result came from the dense fallback. fallback_reason : the repr of the exception that triggered the fallback, or ``None`` when no fallback happened. + residual_fallback : ``True`` when a symmetric path returned modes whose + backward error exceeded + :attr:`~pybmodes.options.SolverOptions.residual_retry_threshold` + and the result was redone through the general dense path. The + symmetric routines factorise the mass matrix, which a very light + beam carrying a very heavy lump makes nearly singular; there they + return wrong low modes rather than failing, so the residual is + what catches it. max_residual : the largest per-mode relative residual ``||K x - λ M x|| / ||K x||`` over the returned modes (``0.0`` when no modes were returned). A healthy modal solve sits near machine precision; a large value flags an ill-conditioned or defective eigenproblem. - residuals : the per-mode relative residuals, one per returned mode. + + Measured against the matrices the returned modes **actually + solve**: the symmetrised pair on the symmetric paths, which + symmetrise internally, and the raw pair on the general one. + Measuring a symmetric solve against the raw matrices would charge + it for the skew it was told to discard — on a model with a wide + dynamic range that reads as a large backward error for a solve + that is exact, which is the opposite of what this field is for. + residuals : the per-mode relative residuals, one per returned mode, + on the same basis as ``max_residual``. matrix_cond : 2-norm condition number of the (symmetrised) mass matrix, or ``None`` when not computed (sparse path, or a system larger than the dense-conditioning size limit). @@ -114,6 +219,7 @@ class SolverDiagnostics: max_residual: float residuals: tuple[float, ...] matrix_cond: float | None + residual_fallback: bool = False # Sparse path activates once the reduced system has more than this # many DOFs and the caller asked for a small subset of modes. Below @@ -223,30 +329,124 @@ def solve_modes( _normalize_columns_l2(eigvecs) + # The residual retry — see the module docstring for the rule and for + # why each of its clauses exists. In short: dense ``eigh`` reduces + # through a Cholesky factor of the *mass* matrix and fails silently + # when that is nearly singular, and the backward error is what + # catches it. + # + # Everything below measures against the matrices the returned modes + # actually solve — the symmetrised pair on a symmetric path, since + # both symmetrise internally. That is passed as a flag rather than by + # building the pair here: materialising it costs two dense ngd-square + # allocations, and a large sparse solve would pay for them on every + # call without ever needing them. + residual_fallback = False + # Dense symmetric only — ``eigsh`` factorises ``K``, so it does not + # have this failure, and its mode window is a different set that must + # not be index-compared. Size-capped because a sparse solve that + # fails to converge falls back to dense at *any* size. Both reasons + # in full in the module docstring. + if ( + sym + and path == "dense_symmetric" + and ngd <= _SOLVER_OPTIONS.residual_retry_max_ndof + ): + # Measured against the matrices the symmetric paths actually + # solved. Both symmetrise internally, and the accepted skew is + # only guaranteed small relative to ``max|K|``: in a model with a + # wide dynamic range it can still be large relative to a soft + # mode's own eigenvalue. Judging an exact symmetric solve against + # the unsymmetrised matrices would then show a residual above the + # threshold, and ``eig`` on those same unsymmetrised matrices + # would "win decisively" purely by answering a different question + # — replacing a correct spectrum with the skew's. + # + # 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) + if sym_r.size and float(sym_r.max()) > _SOLVER_OPTIONS.residual_retry_threshold: + # Now the pair is worth building: ``eig`` needs matrices + # rather than products. Bounded by the size cap above, and + # this branch is rare. + gk_s = 0.5 * (gk + gk.T) + gm_s = 0.5 * (gm + gm.T) + try: + alt_vals, alt_vecs, ordering_sound = _general_spectrum_for_retry( + gk_s, gm_s, n_modes, + ) + except (np.linalg.LinAlgError, ValueError) as exc: + # The alternative is a best-effort second opinion, not a + # requirement. A pencil defective enough to break the + # symmetric reduction can also break ``eig``, and turning + # that into a hard failure would make this guard destroy + # usable results on exactly the inputs it was added to + # help. Decline and keep what we have. + _log.warning( + "solve_modes: residual retry failed (%r); keeping the " + "symmetric result", exc, + ) + alt_vals = np.empty(0) + alt_vecs = np.empty((eigvecs.shape[0], 0)) + ordering_sound = False + if alt_vecs.size: + _normalize_columns_l2(alt_vecs) + alt_r = _modal_residuals(gk_s, gm_s, alt_vals, alt_vecs) + improved, regressed = ( + _compare_candidate_modes( + sym_r, alt_r, alt_vals.size, eigvals.size, + ) + if ordering_sound + else (np.zeros(0, dtype=bool), np.zeros(0, dtype=bool)) + ) + # Accepting replaces the whole spectrum, not just the modes + # that prompted the retry, so a candidate that fixes one mode + # while ruining another is not an improvement to the result. + if improved.any() and not regressed.any(): + idx = int(np.argmax(np.where(improved, sym_r[:improved.size], 0.0))) + warnings.warn( + f"the symmetric eigensolver returned " + f"{int(improved.sum())} mode(s) that do not satisfy " + f"K x = lambda M x — worst at index {idx}, backward " + f"error {sym_r[idx]:.2e} against {alt_r[idx]:.2e} from " + f"the general dense path. The returned modes come from " + f"the general solve, which factorises neither matrix. " + + _retry_cause(gm_s), + RuntimeWarning, + stacklevel=2, + ) + eigvals, eigvecs = alt_vals, alt_vecs + path = "dense_general" + residual_fallback = True + # Mode-count guarantee: the general path filters complex / non- # positive eigenvalues, so it can return fewer modes than requested. # Surface that rather than letting it pass silently (a downstream # broadcast would otherwise fail with an opaque shape error). # - # Gate the warning to the general path only (Codex P2). The dense - # symmetric path also returns fewer than ``n_modes`` when the request - # simply exceeds the available DOFs (it truncates to - # ``min(n_modes, ngd)``), which is a benign "asked for more modes than - # the system has" case, not a defective eigenproblem — warning there - # would mislead, and would fail callers that treat warnings as errors. + # Gate on modes actually *discarded*, not on the path label. Two + # benign shortfalls would otherwise be reported as a defective + # eigenproblem. Asking for more modes than the system has is one: + # every path truncates to ``min(n_modes, ngd)``, which is a request + # the caller can reasonably make. A residual retry is the other — it + # relabels the path ``"dense_general"`` while preserving the whole + # spectrum, so nothing was filtered, and a 117-DOF system asked for + # 1000 modes would be reported as defective for returning its 117. n_returned = int(eigvecs.shape[1]) + n_available = ngd if n_modes is None else min(n_modes, ngd) if ( path == "dense_general" - and n_modes is not None - and n_returned < n_modes + and not residual_fallback + and n_returned < n_available ): warnings.warn( - f"solve_modes recovered only {n_returned} of the requested " - f"{n_modes} modes via the general (non-symmetric) eig path. " - f"The eigenproblem is likely near-degenerate or defective (a " - f"non-symmetric PlatformSupport block can do this); the " - f"missing modes had complex or non-positive eigenvalues and " - f"were filtered out.", + f"solve_modes recovered only {n_returned} of the " + f"{n_available} modes available via the general " + f"(non-symmetric) eig path. The eigenproblem is likely " + f"near-degenerate or defective (a non-symmetric " + f"PlatformSupport block can do this); the missing modes had " + f"complex or non-positive eigenvalues and were filtered out.", RuntimeWarning, stacklevel=2, ) @@ -258,6 +458,7 @@ def solve_modes( gk, gm, eigvals, eigvecs, path=path, symmetric=sym, n_requested=n_modes, sparse_fallback=sparse_fallback, fallback_reason=fallback_reason, + residual_fallback=residual_fallback, ) return eigvals, eigvecs, diagnostics @@ -273,9 +474,19 @@ def _build_diagnostics( n_requested: int | None, sparse_fallback: bool, fallback_reason: str | None, + residual_fallback: bool = False, ) -> SolverDiagnostics: - """Assemble a :class:`SolverDiagnostics` for a completed solve.""" - residuals = _modal_residuals(gk, gm, eigvals, eigvecs) + """Assemble a :class:`SolverDiagnostics` for a completed solve. + + ``symmetric`` selects the basis the residuals are measured on: a + symmetric path solved the symmetrised pair, so charging its modes for + the skew it was told to discard would report a correct solve as + defective. After a residual retry the modes came from ``eig`` on that + same symmetrised pair, so the basis is unchanged. + """ + residuals = _modal_residuals( + gk, gm, eigvals, eigvecs, symmetrise=symmetric, + ) cond = _mass_matrix_cond(gm, path) return SolverDiagnostics( path=path, @@ -287,23 +498,264 @@ def _build_diagnostics( max_residual=float(residuals.max()) if residuals.size else 0.0, residuals=tuple(float(r) for r in residuals), matrix_cond=cond, + residual_fallback=residual_fallback, + ) + + +# Above this the mass matrix is ill-conditioned enough for the Cholesky +# reduction to be the credible culprit; below it, something else in the +# pencil is. +_MASS_COND_ATTRIBUTION = 1.0e8 + + +def _retry_cause(gm: np.ndarray) -> str: + """The explanatory half of the retry warning, attributed honestly. + + The near-singular mass matrix is the *motivating* case, not the only + one: the symmetric reduction degrades on an ill-conditioned pencil + generally, and a stiffness spectrum spanning 1e-16 to 1 with ``M = I`` + triggers this guard while ``cond(M) = 1``. Naming the mass matrix + there would send the reader to check a mass distribution that is + perfectly fine. + + So the cause is measured before it is asserted. The condition number + is only computed on this branch, which is rare, and is skipped for a + system large enough for the O(n^3) estimate to matter — where the + text falls back to naming both possibilities. + """ + if gm.shape[0] > _COND_DENSE_MAX: + return ( + "This happens when the pencil is ill-conditioned — most often " + "a nearly singular mass matrix, from a very light beam " + "carrying a very heavy lump, but a very wide stiffness range " + "does it too. Worth checking the section properties for an " + "extreme mass or stiffness ratio." + ) + try: + cond = float(np.linalg.cond(gm)) + except np.linalg.LinAlgError: + cond = float("inf") + if cond > _MASS_COND_ATTRIBUTION: + return ( + f"The symmetric reduction goes through a Cholesky factor of " + f"the mass matrix, which is nearly singular here " + f"(cond = {cond:.1e}) — a very light beam carrying a very " + f"heavy lump does this. Worth checking the mass distribution " + f"is the one you intended." + ) + return ( + f"The mass matrix is well conditioned (cond = {cond:.1e}), so the " + f"reduction was defeated by the pencil rather than by the mass: a " + f"stiffness range wide enough to put a soft mode at the level of " + f"roundoff will do it. Worth checking the section properties for " + f"an extreme stiffness ratio." ) +def _compare_candidate_modes( + sym_r: np.ndarray, + alt_r: np.ndarray, + n_alt: int, + n_sym: int, +) -> tuple[np.ndarray, np.ndarray]: + """Per-mode verdicts on the alternative: ``(improved, regressed)``. + + Both are needed because accepting the retry replaces the **whole** + spectrum, not the modes that prompted it. A candidate that fixes one + mode while ruining another is not an improvement to the result even + though it is an improvement to that mode, so "some mode got + decisively better" is only half the test; the other half is that no + mode got decisively worse. + + A mode has regressed when it was acceptable, comes back worse, and + lands above the regression floor. Any worsening counts at any ratio: + the improvement side's margin answers a different question — rescue + or noise — and borrowing it here left a mode free to slide from 0.02 + to 0.099 unflagged, which is the edge of tolerance. + + The landing bound is what keeps this usable rather than paralysing: + a mode going from 5.7e-6 to 3.5e-5 is six times worse and three + orders below anything that matters, and flagging it would block + nearly every legitimate rescue. + + A mode already failing in the symmetric solve gets **no regression + verdict**. Above the threshold neither candidate is trustworthy, and + judging that region would let noise veto every rescue that happens to + sit beside a free-free mode. ``max_residual`` still reports it. + + Rigid-body modes are kept out of the *acceptance* side by the + resolution bar rather than by being identified, which three attempts + established cannot be done reliably here — not by eigenvalue scale, + not by strain, and not by which side of the failure threshold the + residual happens to fall on. Dividing one roundoff quantity by + another produces an arbitrary number: 12.4, 0.79 and 0.076 have all + been measured on healthy models, and the last is *below* the failure + threshold, so it looked exactly like a mode being resolved. The one + thing roundoff reliably does not do is land near machine precision. + + The cost is a real case declined: a breakdown the alternative + improves substantially without resolving is not acted on. Neither + result is trustworthy there, so keeping the original and reporting + the backward error is the honest outcome. + + Together the two verdicts give the guarantee the caller relies on: a + mode that was acceptable can only end up above the regression floor + by having *improved*, never as collateral. + + The comparison has to be **per mode**, not on the two maxima. A + rigid-body mode's backward error is a ratio of two near-zero + quantities and reads ~1 in *both* candidates however exact each is, + so it sets a floor under the alternative's maximum: with one present, + ``max(alt_r)`` stays near 1, and no improvement elsewhere can drive + it below the required fraction of ``max(sym_r)`` unless the symmetric + solve is worse still by that same fraction inverted. A free-free + model with a genuinely corrupted elastic mode at a backward error of + ~0.8 would sail through, which is exactly the breakdown this guard + exists to catch. + + Comparing mode by mode removes the floor: the rigid modes contribute + ~1 against ~1 and register as no improvement, while a corrupted + elastic mode contributes ~0.8 against ~1e-9 and registers clearly. + Both candidates are sorted ascending over the same spectrum (the + retry preserves rigid-body modes for this reason), so equal indices + describe the same mode. + + Returns two boolean masks over the compared modes, both empty when + the alternative recovered fewer modes than the symmetric solve — + losing a mode is never an improvement, whatever the residuals say. + """ + empty = np.zeros(0, dtype=bool) + if n_alt < n_sym: + return empty, empty + n = min(sym_r.size, alt_r.size) + if n == 0: + return empty, empty + threshold = _SOLVER_OPTIONS.residual_retry_threshold + factor = _SOLVER_OPTIONS.residual_retry_improvement + sym, alt = sym_r[:n], alt_r[:n] + # What separates a rescue from noise is the *size* of the win, not + # which side of a line the candidate lands on. A rigid-body mode's + # residual divides one near-zero quantity by another, so its value is + # arbitrary — 12.4, 0.79 and 0.076 have all been measured on healthy + # models, and the last is below the failure threshold, so no absolute + # threshold can exclude it. Its *ratio*, though, stays around 11x to + # 16x, while a genuine rescue improves by 1e5 to 1e10. Four orders + # separate the two populations. + # + # The absolute bar is kept as a second condition for the case the + # ratio cannot see: a wildly broken 1e6 against a candidate at 100 + # clears any ratio while both remain garbage. + resolved = _SOLVER_OPTIONS.residual_retry_resolved + improved = (sym > threshold) & (alt <= resolved) & (alt < factor * sym) + # A mode that was acceptable must not come back materially worse. + # Any worsening counts, at any ratio: the improvement side's margin + # answers "is this a rescue or noise", which is a different question, + # and borrowing it here left a mode free to slide from 0.02 to 0.099 + # unflagged. What bounds this instead is where the mode *lands* — + # below the regression floor the change cannot matter, which is what + # keeps harmless churn (5.7e-6 to 3.5e-5) from blocking every rescue. + # + # Modes already failing in the symmetric solve get no verdict at all. + # Neither value is trustworthy there, and a rigid-body mode — whose + # residual divides roundoff by roundoff and has been seen to read + # 12.4 against 0.79 on a healthy pencil — lives entirely in that + # region. Judging it would be judging noise, and doing so in this + # direction would let that noise veto every legitimate rescue. + # ``max_residual`` still reports such a mode to the caller. + floor = _SOLVER_OPTIONS.residual_regression_floor + regressed = (sym <= threshold) & (alt > floor) & (alt > sym) + return improved, regressed + + +# The two routes to ``sym(A) v`` peak at ``3 n k`` and ``2 n^2`` bytes of +# temporaries, so they cross over at ``k = 2 n / 3`` — measured, not +# derived: a first estimate of ``n / 3`` was wrong by a factor of two and +# would have taken the more expensive route across a third of the range. +# Columns per pass of the residual sweep. Every temporary in the sweep +# is ``ngd x`` this, so the peak is bounded by it rather than by the +# number of modes the caller asked for. +# +# Only the peak depends on it — the result does not, since the residual +# is per mode and blocking merely partitions the columns. Chosen at the +# knee of the measured time curve: sweeping the full spectrum of a +# 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 + + +def _prefer_materialised(n_rows: int, n_cols: int) -> bool: + """Is ``sym(A) v`` cheaper built than split into two products? + + The split form peaks at three ``n x k`` temporaries and the built one + at an ``n x n`` copy plus the ``n x k`` result, so they cross over at + ``3 n k = 2 n^2``, i.e. ``k = 2n/3``. + + Measured, not derived. A flop-count estimate put the crossover at + ``n/3``, which would have taken the dearer route across a third of + the range. + """ + return n_cols * 3 >= n_rows * 2 + + +def _apply(a: np.ndarray, v: np.ndarray, symmetrise: bool) -> np.ndarray: + """``A v``, or ``sym(A) v`` by whichever route is cheaper. + + ``0.5 (A + A.T) v == 0.5 (A v + A.T v)``, and the right-hand side + avoids a dense ngd-square allocation — but only while ``v`` is + narrow. At full width the two products cost more than the copy they + were avoiding, so the route is tested rather than assumed. + + Blocking the caller's sweep does not remove the need for the test. + It bounds the block at :data:`_RESIDUAL_BLOCK` columns, which is + narrow relative to a large ``ngd`` but not to a small one: below + ``ngd = 192`` a full-spectrum request still hands this a block wider + than the crossover. Measured at ``ngd = 100``, testing the width + there runs the sweep in 116 us against 178 us for the same peak. + """ + if not symmetrise: + return np.asarray(a @ v) + 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)) + + def _modal_residuals( gk: np.ndarray, gm: np.ndarray, eigvals: np.ndarray, eigvecs: np.ndarray, + *, symmetrise: bool = False, ) -> np.ndarray: """Per-mode relative backward error ``||K x - λ M x|| / ||K x||``. - The honest health metric for a generalised modal solve. Cheap - (matrix-times-thin-matrix), so computed for every path. + The honest health metric for a generalised modal solve, and cheap + enough to compute on every path. + + ``symmetrise`` measures against ``sym(A)`` instead, which is what a + symmetric path actually solved. + + **Swept in column blocks.** The residual is per mode, so no step + needs every mode present at once, and holding them all would tie the + peak to a number the caller chooses: ``n_modes=None`` is the public + default and returns the whole spectrum, which at the retry size cap + would put four ngd-square arrays live at once — over 100 MB — purely + to measure. Blocking bounds every temporary at ``ngd`` by + :data:`_RESIDUAL_BLOCK` instead, and leaves the common thin case + (a handful of modes out of a few thousand DOFs) in a single pass, + computed exactly as before. """ if eigvecs.size == 0: return np.empty(0, dtype=float) - kx = gk @ eigvecs # (ngd, k) - mx = gm @ eigvecs - num = np.linalg.norm(kx - mx * eigvals[np.newaxis, :], axis=0) - den = np.linalg.norm(kx, axis=0) + n_modes_out = eigvecs.shape[1] + num = np.empty(n_modes_out, dtype=float) + den = np.empty(n_modes_out, dtype=float) + for lo in range(0, n_modes_out, _RESIDUAL_BLOCK): + hi = min(lo + _RESIDUAL_BLOCK, n_modes_out) + block = eigvecs[:, lo:hi] + kx = _apply(gk, block, symmetrise) # (ngd, <= block) + mx = _apply(gm, block, symmetrise) + den[lo:hi] = np.linalg.norm(kx, axis=0) + num[lo:hi] = np.linalg.norm( + kx - mx * eigvals[np.newaxis, lo:hi], axis=0, + ) return np.asarray(num / np.where(den > 0.0, den, 1.0), dtype=float) @@ -402,6 +854,57 @@ def _solve_dense_general( return eigvals[order], eigvecs[:, order] +def _general_spectrum_for_retry( + gk: np.ndarray, gm: np.ndarray, n_modes: int | None, +) -> tuple[np.ndarray, np.ndarray, bool]: + """The general solve as the retry in :func:`solve_modes` needs it. + + Separate from :func:`_solve_dense_general` so the asymmetric + production path keeps the BModes-matching filter it is validated + against, untouched. + + Two differences, both about making per-index comparison against a + symmetric solve meaningful. + + **No sign filter.** ``eigh`` filters nothing, so discarding + non-positive eigenvalues here would return a *different set* — the + same length, since a truncated request backfills the gap from higher + up — and equal indices would stop meaning equal modes. Both omissions + are reachable: a free-free model's zero-frequency modes, and the + negative eigenvalues an indefinite ``K`` produces once + ``run(gravity=...)`` loads a column past its buckling weight. + + **A verified ordering.** Complex eigenvalues cannot be kept, and on a + symmetric problem ``eig`` does emit a few rounding-induced conjugate + pairs — routinely, and harmlessly, because they land at the stiff end + of the spectrum far above any mode a caller asks for. Demanding that + none appear is therefore too strict to be useful. What actually + matters is narrower: whether anything discarded would have fallen + *inside* the returned window. The third return value reports that, + and the caller declines to swap when it is ``False``, since an + unverifiable ordering is not a basis for replacing a result. + """ + vals_all, vecs_all = eig(gk, gm) + closed = np.real_if_close(vals_all, tol=1000) + keep_mask = np.isreal(closed) & np.isfinite(closed.real) + + vals = closed.real[keep_mask] + vecs = np.real_if_close(vecs_all[:, keep_mask], tol=1000).real + order = np.argsort(vals) + vals, vecs = vals[order], vecs[:, order] + + keep = vals.size if n_modes is None else min(n_modes, vals.size) + dropped = vals_all[~keep_mask] + ordering_sound = True + if dropped.size and keep: + # Sound exactly when every discarded eigenvalue sits above the + # window, so the window really is the smallest ``keep`` modes. + ordering_sound = bool( + np.min(dropped.real) > vals[keep - 1] + ) + return vals[:keep], vecs[:, :keep], ordering_sound + + # --------------------------------------------------------------------------- # Utility helpers # --------------------------------------------------------------------------- diff --git a/src/pybmodes/options.py b/src/pybmodes/options.py index 075d0e4..09ddfff 100644 --- a/src/pybmodes/options.py +++ b/src/pybmodes/options.py @@ -70,10 +70,74 @@ class SolverOptions: :func:`scipy.linalg.eig` instead of the symmetric :func:`scipy.linalg.eigh`. The OC3 Hywind cross-coupled ``hydro_K + mooring_K`` exercises this branch. + residual_retry_threshold : float, default 0.1 + Largest per-mode relative residual + ``||K x - λ M x|| / ||K x||`` a symmetric solve may return before + the general dense path is tried as well. ``scipy.linalg.eigh`` + reduces the generalised problem through a Cholesky factor of the + **mass** matrix, which loses accuracy once that matrix is nearly + singular — a very light beam carrying a very heavy lump — and + returns confidently wrong low modes rather than failing. + + The band this sits in is narrower than it looks. Ordinary solves + land at or below ~1e-3; a real deck whose adapter leaves ``M`` + genuinely ill-conditioned (the bundled NREL 5MW land tower, cond + ~4e10) reaches ~2e-2 and is *not* meant to trigger; the degraded + regime starts around 0.7. The default splits the last two gaps + with roughly 5x either side. + residual_retry_resolved : float, default 1e-3 + Backward error the candidate must reach on a mode before that + mode can justify a swap. Guards the case the ratio alone cannot: + a wildly broken symmetric solve at 1e6 against a candidate at + 100 clears any ratio while both remain garbage. + residual_regression_floor : float, default 1e-2 + Where a worsened mode has to land before the worsening counts. + Below this the change cannot matter, which is what stops harmless + churn — 5.7e-6 to 3.5e-5 — from blocking every rescue. Kept + separate from ``residual_retry_improvement`` on purpose: they + answer different questions, and deriving one from the other + coupled two unrelated decisions. + residual_retry_max_ndof : int, default 2000 + Largest reduced system the retry will attempt. The retry is a + dense ``eig``, whose cost grows as ``ngd^3`` with a much larger + constant than the ``eigh`` it is checking. Normally that is + bounded by ``sparse_ndof_threshold``, since anything bigger takes + the sparse path and is not retried — but a sparse solve that + *fails to converge* falls back to dense at any size, and there an + unbounded retry could spend minutes on a model whose result was + already available. A guard against a silent wrong answer should + not be able to turn one into a silent hang, so above this size it + declines and leaves the backward error to the diagnostics. + residual_retry_improvement : float, default 1e-3 + How much better the general path's backward error must be on a + mode before that mode can justify taking its result. The + load-bearing guard, and the one that separates a rescue from + noise. + + Two populations were measured while building this. Genuine + rescues improve by 1e5 to 1e10: the worst observed went from + 3.98e1 to 3.17e-4. Rigid-body roundoff, which divides one + near-zero quantity by another and so produces an arbitrary + number, improves by 11x to 16x: 0.848 to 0.0762 on a healthy + model, which no absolute threshold excludes because 0.0762 sits + below the failure line. Four orders separate the two, and the + default sits in the middle with roughly 60x margin on the noise + side and 100x on the rescue side. + + A marginal win is refused for a second reason as well: on the + bundled NREL 5MW land deck the general path is 1.4x better while + *breaking* a physically real degenerate fore-aft / side-side pair + the symmetric solver resolves exactly, which the downstream FA / + SS classifier depends on. """ sparse_ndof_threshold: int = 500 symmetry_rtol: float = 1.0e-12 + residual_retry_threshold: float = 0.1 + residual_retry_improvement: float = 1.0e-3 + residual_retry_resolved: float = 1.0e-3 + residual_regression_floor: float = 1.0e-2 + residual_retry_max_ndof: int = 2000 @dataclass(frozen=True) diff --git a/tests/fem/test_gravity_and_point_mass.py b/tests/fem/test_gravity_and_point_mass.py index 20aae85..4f4f5fe 100644 --- a/tests/fem/test_gravity_and_point_mass.py +++ b/tests/fem/test_gravity_and_point_mass.py @@ -277,13 +277,14 @@ def test_lump_matches_static_cantilever_stiffness(self, tmp_path): a = 60.0 m_lump = 5.0e5 # A beam mass this far below the lump makes the generalised mass - # matrix badly conditioned, which the dense LAPACK subset path - # does not handle well; 101 elements puts the solve on the sparse - # shift-invert path, where it is clean. Nothing about the lump - # placement needs the fine mesh. - tower = _synthetic_tower(tmp_path, mass_den=1.0e-2, n_elements=101) + # matrix nearly singular, which the symmetric eigensolvers handle + # badly. The solver detects that from the backward error and + # redoes the solve on the general path, so a coarse mesh is fine + # here — see tests/fem/test_ill_conditioned_mass.py. + tower = _synthetic_tower(tmp_path, mass_den=1.0e-2, n_elements=25) tower.add_point_mass(a, m_lump) - f = tower.run(6, check_model=False).frequencies[0] + with pytest.warns(RuntimeWarning, match="backward error"): + f = tower.run(6, check_model=False).frequencies[0] f_ref = np.sqrt(3.0 * EI_PHYS / (m_lump * a**3)) / (2.0 * np.pi) assert f == pytest.approx(f_ref, rel=5.0e-3) diff --git a/tests/fem/test_ill_conditioned_mass.py b/tests/fem/test_ill_conditioned_mass.py new file mode 100644 index 0000000..7bd0b62 --- /dev/null +++ b/tests/fem/test_ill_conditioned_mass.py @@ -0,0 +1,1429 @@ +"""The symmetric eigensolvers degrade silently on a near-singular mass +matrix, and the solver has to notice. + +``scipy.linalg.eigh`` reduces ``K x = lambda M x`` through a Cholesky +factor of the **mass** matrix. When that matrix is nearly singular — a +very light beam carrying a very heavy lump — the reduction loses +accuracy, and LAPACK returns confidently wrong low modes rather than +raising. On the case pinned below the dense symmetric path reported +0.103 Hz against a true 0.0436 Hz, a factor of 2.4, with no error and no +warning. + +The sparse path is exempt and must stay exempt: ``eigsh(sigma=0, +mode='normal')`` factorises ``K`` instead, so a near-singular ``M`` does +not degrade it, and its ``which="LM"`` window selects a different set of +modes that must never be compared against the retry's by index. + +The guard is the backward error ``||K x - lambda M x|| / ||K x||``, and +almost every test here exists because some reading of it turned out to +be wrong. Three things it does *not* establish, each learned the hard +way and each pinned below: + +- **A large error is not evidence of a breakdown.** The bundled NREL 5MW + land deck sits at ~2e-2 because its adapter leaves ``M`` at cond ~4e10; + swapping there churns a validated frequency by 0.84 % and splits a + degenerate fore-aft / side-side pair the symmetric solver resolves + exactly, which the FA / SS classifier depends on. +- **A small error is not evidence of a rescue.** A rigid-body mode's + residual divides one roundoff quantity by another, so its value is + arbitrary: 12.4, 0.79 and 0.076 have all been measured on healthy + models, and the last sits *below* the failure threshold. No absolute + bar can separate that from a solved mode. +- **A better mode does not make a better spectrum.** Accepting replaces + every mode, so a candidate that rescues one while pushing another past + the threshold is a trade, not an improvement. + +What does separate the populations is the *size* of the win. Genuine +rescues improve by 1e5 to 1e10; rigid roundoff by 11x to 16x. + +Analytical reference: a cantilever whose beam mass is negligible next to +a tip lump behaves as a spring-mass oscillator on the static tip +stiffness ``3 EI / L^3``, giving ``f = sqrt(3 EI / (m L^3)) / 2 pi`` +(Blevins 1979, Table 8-1). +""" + +from __future__ import annotations + +import warnings + +import numpy as np +import pytest + +from pybmodes.fem.assembly import assemble +from pybmodes.fem.nondim import RM, ROMG, make_params, nondim_tip_mass +from pybmodes.fem.solver import eigvals_to_hz, solve_modes +from pybmodes.io.bmi import TipMassProps + +L = 100.0 +EI = 1.0e10 +M_TIP = 4.0e5 + +# The lump dominates the beam by ~4e5:1 here, which drives cond(M) past +# 1e12 and is what breaks the Cholesky reduction. +LIGHT = 1.0e-2 +REALISTIC = 1.0e3 + + +def _analytic() -> float: + return float(np.sqrt(3.0 * EI / (M_TIP * L**3)) / (2.0 * np.pi)) + + +def _cantilever_with_tip_lump( + nselt: int, mass_den: float, +) -> tuple[np.ndarray, np.ndarray]: + nd = make_params(radius=L, hub_rad=0.0, rot_rpm=0.0) + eiy = EI / nd.ref4 + eli = 1.0 / nselt + el = np.full(nselt, eli) + xb = np.array([1.0 - (i + 1) * eli for i in range(nselt)]) + tip = nondim_tip_mass( + TipMassProps(mass=M_TIP, cm_offset=0.0, cm_axial=0.0, ixx=0.0, + iyy=0.0, izz=0.0, ixy=0.0, izx=0.0, iyz=0.0), + nd, beam_type=2, id_form=1, hub_conn=1, + ) + gk, gm, _ = assemble( + nselt=nselt, el=el, xb=xb, cfe=np.zeros(nselt), + eiy=np.full(nselt, eiy), eiz=np.full(nselt, eiy), + gj=np.full(nselt, 1.0e3 * eiy), eac=np.full(nselt, 100.0), + rmas=np.full(nselt, mass_den / RM), + skm1=np.full(nselt, 1.0e-5), skm2=np.full(nselt, 1.0e-5), + eg=np.zeros(nselt), ea=np.zeros(nselt), omega2=0.0, + sec_loc=np.array([0.0, 1.0]), str_tw=np.zeros(2), hub_conn=1, + tip_mass=tip, + ) + return gk, gm + + +def _first_frequency(nselt: int, mass_den: float): + gk, gm = _cantilever_with_tip_lump(nselt, mass_den) + with pytest.warns(RuntimeWarning, match="backward error"): + eigvals, _vecs, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + return float(eigvals_to_hz(eigvals, ROMG)[0]), diag + + +class TestIllConditionedMassIsCaught: + """Sizes chosen to straddle the dense / sparse dispatch threshold.""" + + @pytest.mark.parametrize("nselt", [13, 27, 53]) + def test_dense_path_recovers_the_analytic_frequency(self, nselt): + f, diag = _first_frequency(nselt, LIGHT) + assert f == pytest.approx(_analytic(), rel=5.0e-3) + assert diag.residual_fallback is True + assert diag.path == "dense_general" + + @pytest.mark.parametrize("nselt", [13, 27, 53]) + def test_the_symmetric_result_would_have_been_wrong(self, nselt): + """Without the guard the answer is not merely imprecise.""" + from scipy.linalg import eigh + + gk, gm = _cantilever_with_tip_lump(nselt, LIGHT) + raw = eigh( + 0.5 * (gk + gk.T), 0.5 * (gm + gm.T), subset_by_index=(0, 3), + )[0] + f_raw = float(eigvals_to_hz(raw, ROMG)[0]) + assert abs(f_raw - _analytic()) / _analytic() > 0.25 + + def test_backward_error_drops_after_the_retry(self): + _f, diag = _first_frequency(27, LIGHT) + assert diag.max_residual < 1.0e-6 + + +class TestHealthyProblemsAreUntouched: + """The guard must not fire on a mass distribution any real structure + would have — the standing rule for a numerical gate.""" + + @pytest.mark.parametrize("nselt", [13, 27, 53]) + def test_realistic_mass_ratio_stays_on_the_symmetric_path(self, nselt): + import warnings as _w + + gk, gm = _cantilever_with_tip_lump(nselt, REALISTIC) + with _w.catch_warnings(): + _w.simplefilter("error") + eigvals, _vecs, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.symmetric is True + f = float(eigvals_to_hz(eigvals, ROMG)[0]) + # The beam mass now matters, so this sits a little below the + # massless-beam closed form rather than on it. + assert 0.9 * _analytic() < f < _analytic() + + def test_healthy_residuals_keep_real_margin_under_the_threshold(self): + """A 400:1 lump-to-beam ratio is already a demanding model and is + the worst healthy backward error in the suite, at ~1.3e-3. The + threshold sits well above that and an order below the + degraded regime, so neither side is close. + + A false positive would cost time rather than accuracy in any + case: the retry keeps whichever of the two solves has the lower + backward error. + """ + from pybmodes.options import DEFAULT_SOLVER_OPTIONS as opt + + worst = 0.0 + for nselt in (13, 27, 53): + gk, gm = _cantilever_with_tip_lump(nselt, REALISTIC) + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + worst = max(worst, diag.max_residual) + assert worst < 0.05 * opt.residual_retry_threshold + + +class TestMarginalImprovementIsRefused: + """The decisive-win condition, pinned directly. + + A synthetic pair where the symmetric solve is above the threshold but + the general path cannot do materially better must keep the symmetric + result — the behaviour that protects the bundled land deck's + degenerate fore-aft / side-side pair. + """ + + def test_marginal_gain_keeps_the_symmetric_result(self, monkeypatch): + import dataclasses + + import pybmodes.fem.solver as solvermod + + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + + # Force the "general path is barely better" situation by making + # the improvement bar unreachable, leaving the threshold tripped. + monkeypatch.setattr( + solvermod, "_SOLVER_OPTIONS", + dataclasses.replace( + solvermod._SOLVER_OPTIONS, residual_retry_improvement=0.0, + ), + ) + import warnings as _w + + with _w.catch_warnings(): + _w.simplefilter("error") + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.path == "dense_symmetric" + + def test_degenerate_pair_comes_back_as_a_pair(self): + """The property the improvement bar exists to protect: a tower + with EI_FA == EI_SS returns its bending modes as one degenerate + pair rather than two separated modes. + + The tolerance is loose on purpose. How exactly the pair resolves + depends on the LAPACK build — this same model splits it at ~1e-16 + on one and ~3e-4 on another — so pinning it tightly tests the + vendor's BLAS rather than pyBmodes. What matters here is that the + two remain the same mode to engineering precision. + """ + gk, gm = _cantilever_with_tip_lump(27, REALISTIC) + eigvals, _v = solve_modes(gk, gm, n_modes=4) + f = eigvals_to_hz(eigvals, ROMG) + assert f[0] == pytest.approx(f[1], rel=1.0e-2) + + +class TestRigidBodyModesAreNotMistakenForBreakdown: + """A free-free model's zero-frequency modes must not corrupt the result. + + For a rigid-body mode ``K x ~ 0`` and ``lambda ~ 0``, so the relative + residual divides one roundoff quantity by another and its value is + arbitrary. Such modes are deliberately *not* identified — three + attempts to do so failed, which the module docstring of + :mod:`pybmodes.fem.solver` records. They are neutralised instead, by + judging the size of the win and by preserving zero eigenvalues so a + false positive is wasteful rather than destructive. + """ + + def _free_free_with_a_zero_mode(self): + """A symmetric platform with no yaw restoring: one exactly zero + eigenvalue among otherwise well-conditioned elastic modes.""" + n = 12 + rng = np.random.default_rng(7) + a = rng.normal(size=(n, n)) + m = a @ a.T + n * np.eye(n) # SPD, well conditioned + b = rng.normal(size=(n, n - 1)) + k = b @ b.T # rank n-1: one zero mode + return 0.5 * (k + k.T), 0.5 * (m + m.T) + + def test_zero_mode_survives_whichever_path_runs(self): + gk, gm = self._free_free_with_a_zero_mode() + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=6, return_diagnostics=True, + ) + # The zero mode survives rather than being filtered away, and the + # spectrum keeps its full width either way. + assert abs(eigvals[0]) < 1.0e-8 * abs(eigvals).max() + assert eigvals.size == 6 + assert diag.n_returned == 6 + + def test_the_metric_really_does_read_about_one_there(self): + """Why the naive raw-maximum reading is untrustworthy, pinned so + the reasoning above stays anchored to a number.""" + from pybmodes.fem.solver import _modal_residuals + + gk, gm = self._free_free_with_a_zero_mode() + eigvals, eigvecs = solve_modes(gk, gm, n_modes=6) + raw = _modal_residuals(gk, gm, eigvals, eigvecs) + rigid = np.argmin(np.abs(eigvals)) + assert raw[rigid] > 0.1 + # Every elastic mode is exact, so the ~1 is the metric failing, + # not the solve. + elastic = np.ones(raw.size, dtype=bool) + elastic[rigid] = False + assert raw[elastic].max() < 1.0e-8 + + def test_the_retry_preserves_zero_eigenvalues(self): + """The property that makes a false positive harmless: the + alternative solve keeps the rigid modes rather than filtering + them, so both candidates describe the same spectrum.""" + from pybmodes.fem.solver import ( + _general_spectrum_for_retry, + _solve_dense_general, + ) + + gk, gm = self._free_free_with_a_zero_mode() + dropped, _v = _solve_dense_general(gk, gm, 6) + kept, _w, _ok = _general_spectrum_for_retry(gk, gm, 6) + assert abs(kept[0]) < 1.0e-8 * float(np.max(np.abs(kept))) + assert np.min(dropped) > 0.0 + # The default filter loses the zero mode and shifts the rest up. + assert kept[1] == pytest.approx(dropped[0], rel=1.0e-6) + + def _six_rigid_dofs(self): + """A model with six genuinely free rigid-body DOFs above a set of + elastic ones — an unmoored floating platform in the limit.""" + n = 14 + rng = np.random.default_rng(11) + a = rng.normal(size=(n, n)) + m = a @ a.T + n * np.eye(n) + b = rng.normal(size=(n, n - 6)) + k = b @ b.T # rank n-6: six zero modes + return 0.5 * (k + k.T), 0.5 * (m + m.T) + + @pytest.mark.parametrize("n_modes", [1, 3, 6]) + def test_a_rigid_only_subset_keeps_its_modes(self, n_modes): + """The case that broke the eigenvalue-scale classifier. + + Every requested mode is rigid-body, so no reference scale drawn + from the subset can say which of them is meaningful. The retry + may well run; what matters is that it cannot + take modes away, because it now preserves zero eigenvalues and + has to win decisively to be accepted at all. + """ + gk, gm = self._six_rigid_dofs() + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=n_modes, return_diagnostics=True, + ) + assert eigvals.size == n_modes + assert np.all(np.abs(eigvals) < 1.0e-8 * float(np.linalg.norm(gk))) + assert diag.n_returned == n_modes + + def test_a_mixed_subset_is_measured_on_its_elastic_modes(self): + """With elastic modes present the metric is meaningful again and + reports them as exact.""" + gk, gm = self._six_rigid_dofs() + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=10, return_diagnostics=True, + ) + assert eigvals.size == 10 + assert np.max(np.abs(eigvals)) > 1.0e-8 * float(np.linalg.norm(gk)) + # Six rigid modes at the bottom, four exact elastic ones above. + assert np.sum(np.abs(eigvals) < 1.0e-12 * np.max(eigvals)) == 6 + assert max(diag.residuals[6:]) < 1.0e-8 + + +class TestRigidModesCannotMaskAnElasticBreakdown: + """A corrupted elastic mode must be caught even when rigid-body modes + sit alongside it. + + Comparing the two candidates on their *maxima* fails here: a + rigid-body mode reads ~1 in both, so it floors the alternative's + maximum and no improvement among the elastic modes can clear a + decisive-win bar unless the symmetric solve is worse than ~10. A + free-free model with an elastic mode corrupted to a backward error of + ~0.8 would then pass silently — the exact failure this guard exists + to catch, reintroduced by the rigid modes' presence. + + Comparing per mode removes the floor. + """ + + def _rigid_plus_ill_conditioned(self, nselt: int = 27): + """Six free rigid DOFs bolted onto the ill-conditioned cantilever, + then rotated so the two blocks are not separable by inspection.""" + gk, gm = _cantilever_with_tip_lump(nselt, LIGHT) + n = gk.shape[0] + big_k = np.zeros((n + 6, n + 6)) + big_m = np.zeros((n + 6, n + 6)) + big_k[:n, :n] = gk + big_m[:n, :n] = gm + # Zero stiffness, unit mass on the six extra DOFs: rigid-body. + big_m[n:, n:] = np.eye(6) * float(np.trace(gm)) / n + q, _ = np.linalg.qr(np.random.default_rng(3).normal(size=(n + 6, n + 6))) + k_rot = q.T @ big_k @ q + m_rot = q.T @ big_m @ q + return 0.5 * (k_rot + k_rot.T), 0.5 * (m_rot + m_rot.T) + + def test_the_maxima_rule_misses_what_the_per_mode_rule_catches(self): + """The mechanism, pinned on residual vectors directly. + + Stated as arithmetic rather than run through LAPACK on purpose: + whether a given matrix pair happens to exhibit the masking + depends on how badly that build's BLAS degrades, which is not + what this is about. These are the numbers the two rules see — + two rigid-body modes reading ~1 in both candidates, one elastic + mode corrupted to 0.8 and fixed to 1e-9, one mode already exact. + """ + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([1.0, 1.0, 0.8, 1.0e-12]) + alt_r = np.array([1.0, 1.0, 1.0e-9, 1.0e-12]) + + # The rigid modes floor the alternative's maximum at ~1, so a + # maxima rule sees no decisive win and leaves the corruption in. + assert not (alt_r.max() < 0.1 * sym_r.max()) + + # Per mode, the corrupted one is unmissable and the rigid ones + # register as exactly what they are: no improvement either way. + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 4, 4) + assert improved.tolist() == [False, False, True, False] + assert not regressed.any() + + def test_a_shorter_alternative_is_never_an_improvement(self): + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([1.0, 0.8, 1.0e-12]) + alt_r = np.array([1.0e-9, 1.0e-9]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 2, 3) + assert not improved.any() + assert not regressed.any() + + def test_the_result_is_never_made_worse(self): + """The portable guarantee when rigid modes and an ill-conditioned + mass matrix coincide. + + Whether the rescue can *fire* here is build-dependent, and + deliberately so. QZ may represent this problem's theoretically + real zero modes as small complex-conjugate pairs that + ``real_if_close`` will not coerce; where those land in the + spectrum differs between LAPACK builds. On one they sit at the + stiff end, far above anything a caller asks for, and the retry + proceeds. On another they are the zero modes themselves, and + keeping the alternative would mean backfilling the gap with + elastic modes and silently shifting the spectrum. + + So the guard verifies its ordering and declines when it cannot, + and the property asserted here is the one that holds either way: + the requested modes come back, in ascending order, and nothing is + lost. ``TestIllConditionedMassIsCaught`` pins the rescue itself on + the case it was built for, which has no rigid modes and behaves + identically everywhere. + """ + gk, gm = self._rigid_plus_ill_conditioned() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=10, return_diagnostics=True, + ) + assert eigvals.size == 10 + assert diag.n_returned == 10 + assert np.all(np.diff(eigvals) >= -1.0e-12 * max(1.0, abs(eigvals).max())) + + def test_an_unverifiable_ordering_declines_the_swap(self, monkeypatch): + """When the alternative drops a mode from inside the window, its + ordering cannot be trusted and the symmetric result stands.""" + import pybmodes.fem.solver as solvermod + + real_fn = solvermod._general_spectrum_for_retry + + def unsound(gk, gm, n_modes): + vals, vecs, _ok = real_fn(gk, gm, n_modes) + return vals, vecs, False + + monkeypatch.setattr(solvermod, "_general_spectrum_for_retry", unsound) + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + with warnings.catch_warnings(): + warnings.simplefilter("error") + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.path == "dense_symmetric" + + +class TestNegativeEigenvaluesSurviveTheRetry: + """An indefinite ``K`` must not have its unstable modes filtered away. + + ``run(gravity=...)`` past a column's buckling weight drives + eigenvalues genuinely negative, and ``eigh`` returns them as-is. If + the alternative solve filtered them, it would come back the same + *length* — the gap backfilled from higher up — while describing a + shifted spectrum, and the per-index comparison would be reading two + different spectra against each other. + """ + + def _indefinite(self): + """The ill-conditioned cantilever with a genuinely negative mode + bolted on and the pair rotated together.""" + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + n = gk.shape[0] + big_k = np.zeros((n + 2, n + 2)) + big_m = np.zeros((n + 2, n + 2)) + big_k[:n, :n] = gk + big_m[:n, :n] = gm + scale = float(np.trace(gk)) / n + big_k[n, n] = -scale # unstable + big_k[n + 1, n + 1] = scale + big_m[n:, n:] = np.eye(2) * float(np.trace(gm)) / n + q, _ = np.linalg.qr(np.random.default_rng(5).normal(size=(n + 2, n + 2))) + k_rot, m_rot = q.T @ big_k @ q, q.T @ big_m @ q + return 0.5 * (k_rot + k_rot.T), 0.5 * (m_rot + m_rot.T) + + def test_the_default_filter_would_drop_the_unstable_mode(self): + from pybmodes.fem.solver import ( + _general_spectrum_for_retry, + _solve_dense_general, + ) + + gk, gm = self._indefinite() + dropped, _v = _solve_dense_general(gk, gm, 6) + kept, _w, _ok = _general_spectrum_for_retry(gk, gm, 6) + assert kept.min() < 0.0 # the unstable mode is present + assert dropped.min() > 0.0 # and absent from the default + # Same length, shifted spectrum — the trap the index comparison + # would otherwise walk into. + assert kept.size == dropped.size + assert kept[1] == pytest.approx(dropped[0], rel=1.0e-6) + + def test_the_symmetric_solve_reports_it_too(self): + """Both paths must agree that the mode exists, or index matching + is meaningless.""" + from scipy.linalg import eigh + + gk, gm = self._indefinite() + w = eigh(gk, gm, subset_by_index=(0, 5))[0] + assert w.min() < 0.0 + + def test_solve_modes_keeps_the_unstable_mode(self): + gk, gm = self._indefinite() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + eigvals, _v = solve_modes(gk, gm, n_modes=6) + assert eigvals.size == 6 + assert eigvals.min() < 0.0 + + +class TestAcceptedSkewDoesNotTriggerTheRetry: + """A tolerated asymmetry must not be read as a solver failure. + + ``_is_effectively_symmetric`` accepts skew up to ``symmetry_rtol`` + times ``max|K|``, and the symmetric paths then solve the symmetrised + matrices. That skew is only small *relative to the largest* entry: in + a model with a wide dynamic range it can be comparable to a soft + mode's own eigenvalue. Measuring the resulting modes against the + unsymmetrised matrices makes an exact solve look broken, and the + general path then "wins decisively" only because it is answering a + different question — returning the skewed spectrum in place of the + symmetric one the caller was promised. + """ + + def _wide_range_with_accepted_skew(self): + from pybmodes.options import DEFAULT_SOLVER_OPTIONS as opt + + # Eigenvalues spanning 1 down to 1e-12, so max|K| is 1 and the + # tolerated skew is ~1e-12 — the same size as the softest mode. + # The symmetry test compares max|A - A.T|, which is twice the + # off-diagonal skew, so stay under half the tolerance. + d = np.array([1.0, 1.0e-6, 1.0e-12]) + gm = np.eye(3) + s = 0.4 * opt.symmetry_rtol * max(1.0, float(np.max(np.abs(d)))) + # Couple the two softest modes: skew between the stiff ones would + # be negligible against their own scale and prove nothing. + gk = np.diag(d) + np.array([[0.0, 0.0, 0.0], + [0.0, 0.0, s], + [0.0, -s, 0.0]]) + return gk, gm, d + + def test_the_pair_is_accepted_as_symmetric(self): + from pybmodes.fem.solver import _is_effectively_symmetric + + gk, gm, _d = self._wide_range_with_accepted_skew() + assert _is_effectively_symmetric(gk) + assert _is_effectively_symmetric(gm) + + def test_no_retry_and_the_symmetric_spectrum_is_returned(self): + gk, gm, d = self._wide_range_with_accepted_skew() + with warnings.catch_warnings(): + warnings.simplefilter("error") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=3, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.path == "dense_symmetric" + # The symmetrised problem's spectrum, not the skewed one. + assert np.allclose(np.sort(eigvals), np.sort(d), rtol=1.0e-6) + + def test_the_diagnostics_report_the_solved_problem_too(self): + """Telemetry meant to be auditable must not charge a correct + solve for the skew it was told to discard.""" + gk, gm, _d = self._wide_range_with_accepted_skew() + with warnings.catch_warnings(): + warnings.simplefilter("error") + _v, _x, diag = solve_modes(gk, gm, n_modes=3, + return_diagnostics=True) + assert diag.residual_fallback is False + assert diag.max_residual < 1.0e-8 + assert max(diag.residuals) < 1.0e-8 + + def test_measuring_against_the_unsymmetrised_pair_would_have_tripped(self): + """The mechanism: the same exact modes look broken when judged + against matrices they were never solved on.""" + from pybmodes.fem.solver import _modal_residuals + + gk, gm, _d = self._wide_range_with_accepted_skew() + gk_s = 0.5 * (gk + gk.T) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + eigvals, eigvecs = solve_modes(gk, gm, n_modes=3) + against_solved = _modal_residuals(gk_s, gm, eigvals, eigvecs) + against_raw = _modal_residuals(gk, gm, eigvals, eigvecs) + assert against_solved.max() < 1.0e-8 + assert against_raw.max() > 0.1 + + +class TestTheRetryVerifiesItsOrderingRatherThanAssumingIt: + """Index matching is only sound if the alternative recovered every + mode. + + A symmetric problem has ``ngd`` real eigenvalues, so if ``eig`` + returns fewer after its real / finite filter, one was discarded — and + a truncated request would have backfilled the gap from higher up, + leaving equal indices pointing at different modes. The sign filter was + one route to that and is switched off; a rounding-induced complex pair + is another that no flag can prevent, so completeness is checked rather + than assumed. + """ + + def test_a_sound_ordering_is_the_precondition_for_swapping(self): + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + with pytest.warns(RuntimeWarning): + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + # This case does swap, so the precondition held. + assert diag.residual_fallback is True + + def test_a_drop_above_the_window_leaves_the_ordering_sound(self): + """Complex pairs at the stiff end are routine and harmless — they + sit far above anything the caller asked for, so the window is + still the smallest ``n_modes``. + + Built from an explicit rotation block rather than by hoping a + real matrix produces one, so it says the same thing on every + LAPACK build. ``[[a, -b], [b, a]]`` has eigenvalues ``a ± bi``; + with ``a`` far above the window and ``b`` far too large for + ``real_if_close`` to coerce, the pair is dropped from well above + the modes being compared. + """ + from pybmodes.fem.solver import _general_spectrum_for_retry + + gm = np.eye(4) + gk = np.array([ + [1.0, 0.0, 0.0, 0.0], + [0.0, 2.0, 0.0, 0.0], + [0.0, 0.0, 1.0e6, -1.0e3], + [0.0, 0.0, 1.0e3, 1.0e6], + ]) + vals, _v, sound = _general_spectrum_for_retry(gk, gm, 2) + assert sound is True + assert vals.size == 2 + assert np.allclose(vals, [1.0, 2.0]) + + def test_a_drop_inside_the_window_makes_the_ordering_unsound(self): + """A discarded eigenvalue below the top of the window means the + window is not the smallest ``n_modes`` and cannot be compared by + index.""" + from pybmodes.fem.solver import _general_spectrum_for_retry + + # Two well-separated real modes plus a complex pair placed below + # them, which no coercion will make real. + gm = np.eye(4) + gk = np.array([ + [1.0, 0.0, 0.0, 0.0], + [0.0, 2.0, 0.0, 0.0], + [0.0, 0.0, 0.0, -1.0e-3], + [0.0, 0.0, 1.0e-3, 0.0], + ]) + vals, _v, sound = _general_spectrum_for_retry(gk, gm, 2) + assert sound is False + assert vals.size == 2 + + +class TestTheWarningAttributesTheCauseHonestly: + """The near-singular mass matrix motivated this guard but is not the + only thing that trips it, and naming it unconditionally sends the + reader to check something that may be perfectly fine.""" + + def test_a_singular_mass_names_the_mass_matrix(self): + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + with pytest.warns(RuntimeWarning, match="nearly singular here"): + solve_modes(gk, gm, n_modes=4) + + def test_a_well_conditioned_mass_does_not_blame_it(self): + """``M = I`` with a stiffness spectrum spanning 1e-16 to 1 trips + the guard at ``cond(M) = 1``.""" + rng = np.random.default_rng(474) + q, _r = np.linalg.qr(rng.normal(size=(3, 3))) + gk = q @ np.diag([1.0e-16, 1.0e-8, 1.0]) @ q.T + gk = 0.5 * (gk + gk.T) + gm = np.eye(3) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + _v, _x, diag = solve_modes(gk, gm, n_modes=3, + return_diagnostics=True) + text = " ".join(str(w.message) for w in caught) + if diag.residual_fallback: + assert "well conditioned" in text + assert "nearly singular" not in text + assert "stiffness ratio" in text + + def test_the_message_states_what_was_measured(self): + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + solve_modes(gk, gm, n_modes=4) + text = " ".join(str(w.message) for w in caught) + assert "do not satisfy" in text + assert "cond =" in text + + +class TestTheRetryIsScopedToTheDensePath: + """Only the dense symmetric path is retried, and that is about which + matrix each routine factorises. + + ``eigh`` reduces through a Cholesky factor of the **mass** matrix, + which is the failure this guard exists for. ``eigsh(sigma=0, + mode='normal')`` factorises ``K`` instead and is unaffected — the + mesh sweep that motivated the work returns correct frequencies on + exactly the meshes large enough to take the sparse path. + + It also removes a spectrum mismatch by construction: ``which="LM"`` + selects the modes nearest zero *in magnitude* while the retry selects + the algebraically smallest, and with negative eigenvalues present + those are different sets that must never be compared by index. + """ + + def test_the_sparse_path_is_left_alone(self, monkeypatch): + import pybmodes.fem.solver as solvermod + + # Force the sparse path on a small problem, then make the + # residual look terrible. The retry must still not run. + gk, gm = _cantilever_with_tip_lump(13, LIGHT) + monkeypatch.setattr(solvermod, "_SPARSE_NDOF_THRESHOLD", 1) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.path == "sparse_shift_invert" + assert diag.residual_fallback is False + + def test_the_sparse_path_gets_the_ill_conditioned_case_right(self): + """Why leaving it alone is safe rather than a gap: it factorises + K, so a near-singular M does not degrade it.""" + gk, gm = _cantilever_with_tip_lump(101, LIGHT) + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.path == "sparse_shift_invert" + f = float(eigvals_to_hz(eigvals, ROMG)[0]) + assert f == pytest.approx(_analytic(), rel=5.0e-3) + + +class TestTheRetryCostIsBounded: + """A guard against a silent wrong answer must not be able to turn one + into a silent hang. + + The retry is a dense ``eig``, whose cost grows as ``ngd^3`` with a + much larger constant than the ``eigh`` it checks. Normally that is + bounded by the sparse dispatch threshold, but a sparse solve that + *fails to converge* falls back to dense at any size. + """ + + def test_a_system_above_the_ceiling_is_not_retried(self, monkeypatch): + import dataclasses + + import pybmodes.fem.solver as solvermod + + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + ngd = gk.shape[0] + monkeypatch.setattr( + solvermod, "_SOLVER_OPTIONS", + dataclasses.replace( + solvermod._SOLVER_OPTIONS, residual_retry_max_ndof=ngd - 1, + ), + ) + called = [] + real = solvermod._general_spectrum_for_retry + monkeypatch.setattr( + solvermod, "_general_spectrum_for_retry", + lambda *a, **k: (called.append(1), real(*a, **k))[1], + ) + with warnings.catch_warnings(): + warnings.simplefilter("error") + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert called == [] + assert diag.residual_fallback is False + + def test_at_the_ceiling_it_still_runs(self, monkeypatch): + import dataclasses + + import pybmodes.fem.solver as solvermod + + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + monkeypatch.setattr( + solvermod, "_SOLVER_OPTIONS", + dataclasses.replace( + solvermod._SOLVER_OPTIONS, + residual_retry_max_ndof=gk.shape[0], + ), + ) + with pytest.warns(RuntimeWarning): + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is True + + +class TestARetryFailureIsNotAHardFailure: + """A pencil defective enough to break the symmetric reduction can + also break ``eig``. Turning that into an exception would make the + guard destroy usable results on exactly the inputs it exists for.""" + + def test_a_raising_retry_keeps_the_symmetric_result(self, monkeypatch): + import pybmodes.fem.solver as solvermod + + def boom(gk, gm, n_modes): + raise np.linalg.LinAlgError("did not converge") + + monkeypatch.setattr(solvermod, "_general_spectrum_for_retry", boom) + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.path == "dense_symmetric" + assert eigvals.size == 4 + + def test_a_value_error_is_handled_the_same_way(self, monkeypatch): + import pybmodes.fem.solver as solvermod + + def boom(gk, gm, n_modes): + raise ValueError("array must not contain infs or NaNs") + + monkeypatch.setattr(solvermod, "_general_spectrum_for_retry", boom) + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + + +class TestARetryThatTradesModesIsRefused: + """Accepting replaces the whole spectrum, so a candidate that fixes + one mode while ruining another is not an improvement to the result. + + The failing shape: the general solve is exact on the modes that + prompted the retry but pushes a previously acceptable mode above the + failure threshold. Judging on "any mode improved" would take it and + hand back a spectrum with a new bad mode in place of an old one. + """ + + def test_a_trade_is_not_an_improvement(self): + from pybmodes.fem.solver import _compare_candidate_modes + + # Mode 0 rescued decisively; mode 3 was fine and is now above the + # threshold and an order worse. + sym_r = np.array([0.52, 1.1e-3, 5.7e-6, 0.017]) + alt_r = np.array([1.1e-10, 3.0e-11, 3.5e-5, 0.26]) + + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 4, 4) + assert improved[0] + assert regressed[3] + # The rule the caller applies. + assert not (improved.any() and not regressed.any()) + + def test_a_small_crossing_of_the_threshold_is_still_a_regression(self): + """A mode sliding from just under the threshold to well over it + worsens by less than the decisive factor, so the mirrored test + alone lets it through. Any crossing counts.""" + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([0.52, 0.09]) + alt_r = np.array([1.0e-10, 0.8]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 2, 2) + assert improved[0] + assert regressed[1] + assert not (improved.any() and not regressed.any()) + + def test_an_exact_mode_driven_to_just_under_the_bar_is_a_regression(self): + """Gating the worsening test on the threshold itself would let a + mode go from machine precision to 0.1 — fifteen orders — while + staying nominally acceptable. It is gated a tenth lower.""" + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([0.52, 1.0e-16]) + alt_r = np.array([1.0e-10, 0.099]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 2, 2) + assert improved[0] + assert regressed[1] + + def test_an_acceptable_mode_cannot_be_degraded_as_collateral(self): + """The guarantee, checked over **every** pair rather than only + the tenfold ones (an earlier version of this test swept only + those, which is how it certified a bound the rule did not hold). + + An earlier version of this test swept only ``alt > 10 * sym``, + which is why it certified a bound the rule did not actually hold: + a mode sliding from 0.02 to 0.099 is under fivefold and lands at + the edge of tolerance, and nothing looked at it. + + The claim is narrow and about the acceptable side only. A mode + that was acceptable can end up above the regression floor + solely by having improved, never as collateral of someone else's + rescue. + """ + from pybmodes.fem.solver import _compare_candidate_modes + from pybmodes.options import DEFAULT_SOLVER_OPTIONS as opt + + t = opt.residual_retry_threshold + bound = opt.residual_regression_floor + grid = np.logspace(-16, 2, 37) + for s in grid: + if s > t: + continue # not the acceptable side + for a in grid: + imp, reg = _compare_candidate_modes( + np.array([s]), np.array([a]), 1, 1, + ) + if reg[0] or imp[0]: + continue + assert a <= bound or a <= s, ( + f"sym={s:.2e} -> alt={a:.2e} escaped unflagged above " + f"the {bound:.2e} bound without improving" + ) + + def test_a_mode_already_failing_gets_no_verdict(self): + """The deliberate hole, stated so it is not mistaken for one. + + Above the threshold neither candidate is trustworthy, and a + rigid-body mode lives entirely there. Judging that region either + way turns roundoff into a decision. + """ + from pybmodes.fem.solver import _compare_candidate_modes + from pybmodes.options import DEFAULT_SOLVER_OPTIONS as opt + + t = opt.residual_retry_threshold + for sym, alt in [(0.2, 5.0), (12.39, 0.794), (0.794, 12.39), (1.0, 1.0)]: + assert sym > t + _imp, reg = _compare_candidate_modes( + np.array([sym]), np.array([alt]), 1, 1, + ) + assert not reg[0], f"sym={sym} alt={alt} should carry no verdict" + + def test_improved_and_regressed_are_mutually_exclusive(self): + """A single mode cannot be both, or the caller's rule would be + reading a contradiction.""" + from pybmodes.fem.solver import _compare_candidate_modes + + grid = np.logspace(-16, 2, 37) + for s in grid: + for a in grid: + imp, reg = _compare_candidate_modes( + np.array([s]), np.array([a]), 1, 1, + ) + assert not (imp[0] and reg[0]), f"sym={s:.2e} alt={a:.2e}" + + def test_a_mode_that_worsens_but_stays_acceptable_is_not_a_regression(self): + """Mode 2 above goes from 5.7e-6 to 3.5e-5 — six times worse and + entirely irrelevant, since it is nowhere near the threshold.""" + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([0.52, 5.7e-6]) + alt_r = np.array([1.1e-10, 3.5e-5]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 2, 2) + assert improved[0] + assert not regressed.any() + + def test_rigid_body_noise_is_not_a_regression(self): + """Residuals that read ~1 in both candidates wobble either way. + A bare ``alt > sym`` test would call that a regression and block + every rescue that happens to sit beside a free-free mode.""" + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([1.0, 0.8]) + alt_r = np.array([1.0001, 1.0e-9]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 2, 2) + assert improved[1] + assert not regressed.any() + + @pytest.mark.parametrize("alt_rigid", [12.39, 0.794, 0.0762, 0.05]) + def test_rigid_body_noise_cannot_justify_a_swap(self, alt_rigid): + """Rigid residuals divide roundoff by roundoff, so the value is + arbitrary — 12.39, 0.794 and 0.0762 have all been measured on + healthy models, and the last two sit *below* the failure + threshold. No threshold can exclude them; a resolution bar near + machine precision can, because roundoff does not land there. + """ + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([0.848, 3.0e-15, 2.0e-15]) + alt_r = np.array([alt_rigid, 1.0e-15, 2.0e-15]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 3, 3) + assert not improved.any() + assert not regressed.any() + + def test_the_two_populations_are_separated_by_the_ratio(self): + """The measurement the rule is calibrated on, kept as a test so + the constants cannot drift away from their evidence. + + The mode that justifies a swap is judged by *how much* the + candidate improves it, because a rigid residual's value is + arbitrary while its ratio is not. Genuine rescues improve by 1e5 + to 1e10; rigid roundoff by 11x to 16x. + """ + from pybmodes.options import DEFAULT_SOLVER_OPTIONS as opt + + rescues = [(1.56e0, 4.10e-10), (3.01e0, 1.96e-09), + (7.55e1, 4.93e-09), (3.98e1, 3.17e-04)] + noise = [(12.39, 0.794), (0.848, 0.0762)] + + worst_rescue = max(a / s for s, a in rescues) + best_noise = min(a / s for s, a in noise) + assert worst_rescue < opt.residual_retry_improvement < best_noise + # And the absolute bar admits every rescue. + assert max(a for _s, a in rescues) <= opt.residual_retry_resolved + + @pytest.mark.parametrize("seed", [14946, 7, 101, 2024, 55555]) + def test_healthy_free_free_pencils_never_trigger_a_swap(self, seed): + """The empirical half of the argument, over several draws rather + than one. A rank-deficient K with a well-conditioned M is a + healthy free-free model; whatever its null-mode roundoff happens + to read, it must never replace the spectrum.""" + rng = np.random.default_rng(seed) + n = 3 + a = rng.normal(size=(n, n)) + gm = a @ a.T + n * np.eye(n) + b = rng.normal(size=(n, n - 1)) + gk = b @ b.T + gk, gm = 0.5 * (gk + gk.T), 0.5 * (gm + gm.T) + assert np.linalg.cond(gm) < 100.0 # genuinely healthy + with warnings.catch_warnings(): + warnings.simplefilter("error") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=n, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert eigvals.size == n + + def test_a_healthy_free_free_pencil_is_left_alone(self): + """End to end on the shape from the report: rank-deficient K, a + well-conditioned M, exact elastic modes.""" + rng = np.random.default_rng(19) + a = rng.normal(size=(4, 4)) + gm = a @ a.T + 4.0 * np.eye(4) + b = rng.normal(size=(4, 3)) + gk = b @ b.T + gk, gm = 0.5 * (gk + gk.T), 0.5 * (gm + gm.T) + with warnings.catch_warnings(): + warnings.simplefilter("error") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert eigvals.size == 4 + + def test_an_improvement_that_leaves_the_mode_failing_is_declined(self): + """The cost of the rule, stated rather than hidden: a hundredfold + gain that still ends above the threshold is not acted on, because + neither result is trustworthy there.""" + from pybmodes.fem.solver import _compare_candidate_modes + + sym_r = np.array([50.0]) + alt_r = np.array([0.5]) + improved, regressed = _compare_candidate_modes(sym_r, alt_r, 1, 1) + assert not improved.any() + assert not regressed.any() + + def test_end_to_end_a_trading_candidate_is_declined(self, monkeypatch): + import pybmodes.fem.solver as solvermod + + gk, gm = _cantilever_with_tip_lump(27, LIGHT) + real = solvermod._general_spectrum_for_retry + + def trading(gk_, gm_, n_modes): + vals, vecs, sound = real(gk_, gm_, n_modes) + # Corrupt the last returned mode so it is decisively worse. + vecs = vecs.copy() + vecs[:, -1] = np.roll(vecs[:, -1], 1) + return vals, vecs, sound + + monkeypatch.setattr( + solvermod, "_general_spectrum_for_retry", trading, + ) + with warnings.catch_warnings(): + warnings.simplefilter("error") + _v, _x, diag = solve_modes( + gk, gm, n_modes=4, return_diagnostics=True, + ) + assert diag.residual_fallback is False + assert diag.path == "dense_symmetric" + + +class TestTheModeCountWarningStaysHonest: + """A shortfall is only newsworthy when modes were actually discarded. + + Two benign cases would otherwise be reported as a defective + eigenproblem: asking for more modes than the system has, which every + path truncates, and a residual retry, which relabels the path + ``"dense_general"`` while preserving the whole spectrum. + """ + + def test_an_overlarge_request_after_a_retry_is_not_reported(self): + gk, gm = _cantilever_with_tip_lump(13, LIGHT) + ngd = gk.shape[0] + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + eigvals, _v, diag = solve_modes( + gk, gm, n_modes=ngd + 500, return_diagnostics=True, + ) + assert diag.residual_fallback is True + assert eigvals.size == ngd + assert not any( + "recovered only" in str(w.message) for w in caught + ), [str(w.message) for w in caught] + + def test_an_overlarge_request_without_a_retry_is_not_reported(self): + gk, gm = _cantilever_with_tip_lump(13, REALISTIC) + ngd = gk.shape[0] + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + eigvals, _v = solve_modes(gk, gm, n_modes=ngd + 500) + assert eigvals.size == ngd + assert not any("recovered only" in str(w.message) for w in caught) + + def test_a_genuine_shortfall_is_still_reported(self): + """An asymmetric pencil whose modes really are filtered away must + still say so.""" + # Antisymmetric K: every eigenvalue is imaginary, so the general + # path recovers none of them. + gk = np.array([[0.0, -1.0], [1.0, 0.0]]) + gm = np.eye(2) + with pytest.warns(RuntimeWarning, match="recovered only"): + solve_modes(gk, gm, n_modes=2) + + +def _prefer_materialised_for(n_rows: int, n_cols: int) -> bool: + """Thin wrapper so the width-rule cases read as arithmetic.""" + from pybmodes.fem.solver import _prefer_materialised + + return _prefer_materialised(n_rows, n_cols) + + +class TestTheGuardDoesNotTaxEverySolve: + """Residuals are measured on every eligible solve, healthy ones + included, so measuring them must not cost more than it saves. + + Two things keep it cheap. ``0.5 (A + A.T) v == 0.5 (A v + A.T v)`` + and only the second form avoids a dense ngd-square pair, which + matters most on the large sparse solves the copy would hurt. And the + sweep runs in column blocks, so its peak is set by the block width + rather than by how many modes the caller asked for — without that, + the identity above inverts at full width and ``n_modes=None`` is the + public default. + """ + + def _pair(self, n): + rng = np.random.default_rng(3) + a = rng.normal(size=(n, n)) + gm = a @ a.T + n * np.eye(n) + b = rng.normal(size=(n, n)) + gk = b @ b.T + return 0.5 * (gk + gk.T), 0.5 * (gm + gm.T) + + def test_the_two_mechanisms_are_not_interchangeable(self): + """Blocking and the width test answer different questions, and + neither subsumes the other. Capping the block to the crossover + instead of testing the width was measured and is worse at both + ends: at ngd = 100 it splits one BLAS call into two for 206 us + against 116 us, and at ngd = 1000 it lifts the peak from 5.1 MB + to 21.3 MB.""" + from pybmodes.fem import solver + + # Large ngd: blocking is what bounds the peak, and the block it + # produces is already narrow, so the width test does not fire. + assert not _prefer_materialised_for(2000, solver._RESIDUAL_BLOCK) + # Small ngd: blocking barely bites, and the width test does. + assert _prefer_materialised_for(100, solver._RESIDUAL_BLOCK) + + def test_symmetrised_products_match_materialising(self): + from pybmodes.fem.solver import _modal_residuals + + gk, gm = self._pair(40) + rng = np.random.default_rng(5) + v = np.linalg.qr(rng.normal(size=(40, 4)))[0] + w = np.linspace(1.0, 2.0, 4) + ks, ms = 0.5 * (gk + gk.T), 0.5 * (gm + gm.T) + assert np.allclose( + _modal_residuals(gk, gm, w, v, symmetrise=True), + _modal_residuals(ks, ms, w, v), + rtol=1.0e-12, atol=1.0e-15, + ) + + def test_an_asymmetric_pair_is_measured_unsymmetrised(self): + """The flag must actually change the basis, not be inert.""" + from pybmodes.fem.solver import _modal_residuals + + gk = np.array([[1.0, 0.9], [0.0, 2.0]]) + gm = np.eye(2) + rng = np.random.default_rng(9) + v = np.linalg.qr(rng.normal(size=(2, 2)))[0] + w = np.array([1.0, 2.0]) + raw = _modal_residuals(gk, gm, w, v) + sym = _modal_residuals(gk, gm, w, v, symmetrise=True) + assert not np.allclose(raw, sym) + + def test_the_threshold_check_reads_the_caller_s_own_matrices( + self, monkeypatch, + ): + """The check runs on every eligible solve, healthy ones included, + before the threshold has been tested. It must therefore read the + arrays it was given rather than a symmetrised copy of them. + + Asserted by identity rather than by measuring memory: the solve + legitimately allocates elsewhere — ``eigh`` takes matrices, so + ``_solve_dense_symmetric`` builds the pair it needs — which would + swamp any peak-usage threshold and make the test meaningless. + """ + import pybmodes.fem.solver as solvermod + + seen = [] + real = solvermod._modal_residuals + + def spy(k, m, vals, vecs, *, symmetrise=False): + seen.append((k is gk, m is gm, symmetrise)) + return real(k, m, vals, vecs, symmetrise=symmetrise) + + monkeypatch.setattr(solvermod, "_modal_residuals", spy) + gk, gm = self._pair(60) + _v, _x, diag = solve_modes(gk, gm, n_modes=4, return_diagnostics=True) + + assert diag.path == "dense_symmetric" + assert diag.residual_fallback is False + assert seen, "the threshold check did not run" + # Every call took the caller's arrays and asked for the + # symmetrised basis, rather than being handed a built pair. + assert all(is_gk and is_gm and sym for is_gk, is_gm, sym in seen), seen + + @pytest.mark.parametrize("k_over_n", [0.01, 0.2, 0.5, 1.0]) + def test_no_pass_of_the_sweep_exceeds_the_block( + self, k_over_n, monkeypatch, + ): + """The invariant the memory bound rests on: the peak follows the + block width, not the caller's ``n_modes`` — which defaults to + ``None`` and asks for the whole spectrum. Measured at ngd = 2000 + with every mode requested, the unblocked sweep peaked at 128 MB + against 10.3 MB blocked. + + This bounds the peak. It does *not* make the block narrow + relative to ``ngd``, which is a separate question settled by the + width test inside ``_apply`` — see + :meth:`test_a_small_system_still_gets_the_width_test`. + + Deliberately not a ``tracemalloc`` comparison. That reports what + the allocator did, and a freed block from an earlier call gets + reused by a later one, so competing forms time-share buffers + differently by measurement order and by platform. An earlier + version of this test compared three live peaks and failed on + Linux while passing on Windows, having measured reuse. + """ + from pybmodes.fem import solver + + widths = [] + real = solver._apply + monkeypatch.setattr( + solver, "_apply", + lambda a, v, s: (widths.append(v.shape[1]), real(a, v, s))[1], + ) + + n = 200 + gk, gm = self._pair(n) + k = max(1, round(k_over_n * n)) + rng = np.random.default_rng(12) + v = np.linalg.qr(rng.normal(size=(n, k)))[0] + solver._modal_residuals( + gk, gm, np.linspace(1.0, 2.0, k), v, symmetrise=True, + ) + + assert widths, "the sweep did not run" + assert max(widths) <= solver._RESIDUAL_BLOCK, ( + f"k={k}: a pass saw {max(widths)} columns, above the " + f"{solver._RESIDUAL_BLOCK}-column block" + ) + assert sum(widths) == 2 * k, "every mode is measured exactly once" + + def test_a_thin_request_still_runs_in_one_pass(self, monkeypatch): + """Blocking must not tax the common case — a handful of modes out + of a few thousand DOFs is one pass per matrix, exactly the two + products the unblocked form did.""" + from pybmodes.fem import solver + + calls = [] + real = solver._apply + monkeypatch.setattr( + solver, "_apply", + lambda a, v, s: (calls.append(v.shape[1]), real(a, v, s))[1], + ) + + gk, gm = self._pair(400) + rng = np.random.default_rng(14) + v = np.linalg.qr(rng.normal(size=(400, 6)))[0] + solver._modal_residuals( + gk, gm, np.linspace(1.0, 2.0, 6), v, symmetrise=True, + ) + assert calls == [6, 6], f"expected one pass per matrix, got {calls}" + + @pytest.mark.parametrize("k", [1, 7, 128, 129, 260]) + def test_blocking_does_not_change_the_answer(self, k): + """Straddles the block boundary: one pass, one plus a remainder, + and several.""" + from pybmodes.fem.solver import _apply, _modal_residuals + + n = 300 + gk, gm = self._pair(n) + rng = np.random.default_rng(13) + v = np.linalg.qr(rng.normal(size=(n, k)))[0] + w = np.linspace(1.0, 2.0, k) + + kx, mx = _apply(gk, v, True), _apply(gm, v, True) + den = np.linalg.norm(kx, axis=0) + expected = np.linalg.norm( + kx - mx * w[np.newaxis, :], axis=0, + ) / np.where(den > 0.0, den, 1.0) + assert np.allclose( + _modal_residuals(gk, gm, w, v, symmetrise=True), expected, + rtol=1.0e-12, atol=1.0e-15, + ) + + @pytest.mark.parametrize("n", [30, 60, 100, 150]) + def test_a_small_system_still_gets_the_width_test(self, n, monkeypatch): + """Blocking bounds the block at 128 columns, which is narrow + against a large ``ngd`` but not a small one. Below ``ngd = 192`` + a full-spectrum request hands ``_apply`` a block wider than the + ``2n/3`` crossover, so the width test inside it is still load + bearing — dropping it cost 178 us against 116 us at ngd = 100, + for the same peak. + """ + from pybmodes.fem import solver + + assert n < 192, "the point of this case is that 128 is not narrow" + k = n # the n_modes=None default + assert min(solver._RESIDUAL_BLOCK, k) >= 2 * n / 3, ( + "this ngd no longer produces a wide block; pick a smaller one" + ) + assert _prefer_materialised_for(n, min(solver._RESIDUAL_BLOCK, k)) + + @pytest.mark.parametrize("n_rows,n_cols,built", [ + (400, 400, True), # full spectrum, the wide case + (2000, 2000, True), + (1500, 6, False), # a few modes out of many DOFs + (400, 4, False), + (100, 128, True), # small ngd: the block is not narrow + (400, 128, False), # large ngd: the same block is + ]) + def test_the_route_is_the_cheaper_of_the_two(self, n_rows, n_cols, built): + """Asserted against the cost model — ``3nk`` split against + ``2n^2`` built — rather than a measured peak, for the + platform-independence reason given above.""" + assert _prefer_materialised_for(n_rows, n_cols) is built + split, build = 3 * n_rows * n_cols, 2 * n_rows * n_rows + assert _prefer_materialised_for(n_rows, n_cols) == (build <= split) + + @pytest.mark.parametrize("block", [1, 3, 16, 64, 1000]) + def test_the_block_width_is_a_free_parameter(self, block, monkeypatch): + """It may be tuned for the peak-versus-overhead trade without + touching any answer. The residual is per mode, so blocking only + partitions independent columns — a block width that changed the + result would mean the columns were not independent after all.""" + from pybmodes.fem import solver + + n, k = 120, 37 + gk, gm = self._pair(n) + rng = np.random.default_rng(15) + v = np.linalg.qr(rng.normal(size=(n, k)))[0] + w = np.linspace(1.0, 2.0, k) + + monkeypatch.setattr(solver, "_RESIDUAL_BLOCK", 128) + reference = solver._modal_residuals(gk, gm, w, v, symmetrise=True) + monkeypatch.setattr(solver, "_RESIDUAL_BLOCK", block) + assert np.allclose( + solver._modal_residuals(gk, gm, w, v, symmetrise=True), + reference, rtol=1.0e-12, atol=1.0e-15, + ) + + @pytest.mark.parametrize("k", [1, 7, 26, 27, 40]) + def test_both_routes_agree_numerically(self, k): + """Whichever route is taken, the answer is ``sym(A) v``.""" + from pybmodes.fem.solver import _apply + + n = 40 + rng = np.random.default_rng(6) + a = rng.normal(size=(n, n)) + v = rng.normal(size=(n, k)) + assert np.allclose(_apply(a, v, True), 0.5 * (a + a.T) @ v) + + def test_measuring_does_not_allocate_a_matrix_copy(self): + import tracemalloc + + from pybmodes.fem.solver import _modal_residuals + + n = 400 + gk, gm = self._pair(n) + rng = np.random.default_rng(7) + v = np.linalg.qr(rng.normal(size=(n, 5)))[0] + w = np.linspace(1.0, 2.0, 5) + + tracemalloc.start() + base = tracemalloc.get_traced_memory()[0] + _modal_residuals(gk, gm, w, v, symmetrise=True) + peak = tracemalloc.get_traced_memory()[1] - base + tracemalloc.stop() + + # Well under a single dense copy, which is what materialising the + # symmetrised pair would have cost twice over. + assert peak < 0.25 * gk.nbytes, ( + f"peak {peak / 1e6:.1f} MB against a {gk.nbytes / 1e6:.1f} MB " + f"matrix — the symmetrised pair is being materialised" + ) + + +class TestDiagnosticsContract: + def test_residual_fallback_defaults_to_false(self): + gk, gm = _cantilever_with_tip_lump(13, REALISTIC) + _v, _x, diag = solve_modes(gk, gm, n_modes=4, return_diagnostics=True) + assert diag.residual_fallback is False + + def test_plain_two_tuple_return_still_works(self): + """The retry must not change the historical return shape.""" + gk, gm = _cantilever_with_tip_lump(13, LIGHT) + with pytest.warns(RuntimeWarning): + out = solve_modes(gk, gm, n_modes=4) + assert isinstance(out, tuple) + assert len(out) == 2