evd: scale the input into range as ?GEEV does (#112)

PATCHSET.md

@@ -101,6 +101,7 @@ gate has changed and two kinds of adoption are now in the ledger.
 | Native gemm thread budget in flops | `matmul_imp` handed `private_gemm_x86::gemm` the caller's thread count above `4096 * size_of::<T>()` MACs, and that kernel threads whatever it is handed (unlike the `gemm` crate's kernels, which stop at `48 * 48 * 256`). f64 `32^3` (65 ki flops) was 2.31 us serial and 12.02 us at eight threads; f64 `64^3` kept its best time but lost 1.2-1.4x in the median. The budget is now handed over only above the `gemm` crate's own floor in flops (`2 * gemm::get_threading_threshold()`), counted as flops so complex keeps its two-times gain at `64^3`. After: f64 `32^3` 2.30 us at every budget, `64^3` 17.0 (was 18.3 best / 21.1 median at 8T), `96^3` and c64 unchanged. | KEEP | — | Probe: `docs/worklogs/2026-10-18-small-gemm-thread-floor.md`; campaign pair `20261018T014801Z-baseline-gemm-guard-mkl` / `20261018T015*Z-candidate-gemm-thread-floor-mkl`. Test: `thread_budget_does_not_change_the_product`. |
 | A diagonal product reads the depth in either direction (#113) | `spicy_matmul` with a diagonal `D` (`A D B`) hands its operands to `t4a-private-gemm-x86`, which turns a negative depth stride of its left operand around by reversing that operand's columns and the right operand's rows. Up to 0.1.26 it did not reverse `D`, so each depth index met the wrong `D` entry. The unpivoted LDL's blocked update multiplies `L10 D0 L10^H` this way, so a matrix stored with reversed rows and columns factored wrongly at n >= 65 (LAPACK indefinite family, layouts 2 and 3: ratios of 1e9 to 1e15 in every dtype, positive definite input included). #117 reversed the depth of `A`, `B` and `D` together in `spicy_matmul` as a workaround. The kernel now reverses the diagonal itself (0.1.27, tensor4all/t4a-private-gemm-x86#5, fixing #4 there), so faer requires 0.1.27 and the workaround is removed (#119); the regression test stays. Upstream has the same defect: `private-gemm-x86` 0.1.22, the latest on crates.io (and 0.1.20, the version faer 0.24.4 requires), reverses `lhs` and `rhs` in `src/lib.rs:1244` without `real_diag`. Upstream faer 0.24.4 passes `D` with the operands' raw strides (`src/linalg/matmul/internal/mod.rs:192-194`), so its unpivoted LDL of a matrix stored with reversed rows and columns is affected the same way. sarah-ek/gemm (`gemm-common`, `gemm-f64` 0.19.0) has no product with a diagonal and is not involved. Not reported upstream. | KEEP | tensor4all/faer#113, tensor4all/faer#119 | Test `a_diagonal_product_reads_the_depth_in_any_direction`: both operands in both directions, column- and row-major destination, exact against the plain sum. It fails against 0.1.26 without the workaround, at `B` reversed into a row-major destination, and passes against 0.1.27 without it. The LAPACK `indefinite` family (`lapack_indefinite_*`, `--include-ignored`) passes in all four dtypes. Record: [worklog](docs/worklogs/2026-10-29-diagonal-product-depth-direction.md). Upstream report: [`upstream-reports/gemm-x86-001-diag-negative-stride`](upstream-reports/gemm-x86-001-diag-negative-stride/). |
 | Eigenvector back substitution clamps relative to the eigenvalue (#112) | `evd_from_real_schur_imp` and `evd_from_cplx_schur_imp` replace a vanishing shift `t_ii - lambda` by a floor. The floor was `eps * sum |t_ij|` over the whole Schur form, which grows with `n`. A clamped row leaves a residual of the floor's size in the eigenvector, so the GET22 residual ratio grew like `n` on defective and tightly clustered spectra: `(2n - 1) / 2` on a Jordan block, 62.5 at n = 63, where OpenBLAS `dgeev` gives 0.5. The floor is now `eps * max(|lambda|, max |t_ij|)` (`|re| + |im|`, scaled by `eps` before summing so that it does not overflow), at least the smallest positive subnormal. That is `?TREVC`'s `ulp * (|wr| + |wi|)`, with the entries standing in for `?GEEV`'s prescaling. A clamp this small is range safe only on entries near one, so a Schur form whose largest entry is outside `[2^-20, 2^20]` is first scaled by a power of two, as `?GEEV` scales its input. Without that, a Schur form near the top of the range overflowed in the substitution's products, and one near the bottom had a clamp whose complex reciprocal overflows. The LAPACK `nonsymmetric_evd` diagnostics drop from 306 / 282 / 220 / 216 to 192 / 168 / 112 / 116 (f64 / f32 / c32 / c64). Kinds 3 and 5 pass. On the remaining fixtures that were compared (kinds 2, 4, 6 and 7 at n = 63, kind 4 at 129), LAPACK's own `dgeev` gives 17.5 to 34.4, against CHKHS's threshold of 20, which is calibrated at n <= 16. The comparison does not show a faer-specific defect there, but it does not cover every remaining fixture, side and precision. The threshold is left to the maintainer. Upstream has the same defect: faer 0.24.4 `src/linalg/evd/mod.rs` sums `|a_ij|` over the Schur form (`:567-572` real, `:632-637` complex) and clamps to `eps * norm` (`:701`, `:794`, `:828`, `:917`). On a Jordan block its `Mat::eigen` gives a GET22 ratio of 62.61 at n = 63, against 0.50 on the fork; record and reproduction in [`upstream-reports/faer-003-evd-clamp-grows-with-n`](upstream-reports/faer-003-evd-clamp-grows-with-n/). Not reported upstream. | KEEP | tensor4all/faer#112 | Test `lapack_nonsymmetric_evd_clustered_residuals_do_not_grow_with_n` (kinds 3 and 5, n = 63 and 65, both sides, residual < 20, every dtype) runs with the regular tests and fails without the change (62.5). `jordan_chains_at_both_ends_of_the_range` covers a Jordan chain scaled to 2^972 and 2^-972 (2^105 and 2^-105 for f32 and c32), and fails without the prescaling. The existing range-limit and subnormal eigenvector tests pass. Record: [worklog](docs/worklogs/2026-10-29-eigenvector-clamp.md). Upstream report: [`upstream-reports/faer-003-evd-clamp-grows-with-n`](upstream-reports/faer-003-evd-clamp-grows-with-n/). |
+| The nonsymmetric EVD scales its input as `?GEEV` does (#112) | `evd_imp` ran the Hessenberg reduction and the Schur iteration on the matrix as given. `?GEEV` first scales it into `[smlnum, bignum]`, `smlnum = sqrt(safmin) / eps`. Without that, a matrix with entries near the bottom of the range was reduced and iterated in subnormals: on random, Grcar and near-Jordan matrices the GET22 ratio grew from about `2^-90` (32-bit) and was `1e4` to `2e15` from `2^-100` (32-bit) and `2^-960` (64-bit) down, where LAPACK gives 2 to 70, and near the top the complex Schur iteration kept the `c32` LAPACK group failing (kinds 6 and 7 at n = 63 to 65: 41.7 to 59.4 against `cgeev`'s 29 to 40; Schur backward error 226 to 258, 67 to 93 when scaled). `evd_imp` now scales by a power of two to just inside the bound (no rounding of its own) and scales the eigenvalues back. The bound is tested on the largest `|re|` or `|im|`, because faer's complex modulus overflows for `2^1023 + 0i` and underflows for `2^-1074 + 0i`. | KEEP | tensor4all/faer#112 | Test `geev_scaling_tests::evd_scales_its_input_into_range` (Grcar and pseudo-random n = 32, f32/c32 at `2^-120` to `2^120`, f64/c64 at `2^-1010` to `2^1015`, right vectors of `a` and `a^H`, GET22 below `20 * max(1, sqrt(n / 16))`) fails without the scaling (`f32`, `2^-120`: 1.7e6) and passes with it; `evd_scaling_at_the_complex_range_ends` (`2^1023 + 0i`, `2^127 + 0i`, a rotation by the smallest subnormal) returns the exact eigenvalues and fails with a modulus-based scan (inf modulus, zero factor, NaN eigenvalue with `Ok`). The `c32` nonsymmetric EVD LAPACK group passes, so its known-failures entry is removed; all LAPACK validation families pass. A 372-matrix sweep (f32/f64/c32/c64, n = 16 to 257, random, Grcar, near-Jordan, graded, scales to both range limits) against LAPACK `?geev`: every non-graded series is at most 1.74 times LAPACK after the change. It also found #144 (the real Schur deflation sort can loop forever) and #145 (real eigenvectors of graded matrices). The arithmetic count is unchanged up to one pass over the input. Upstream report: [`upstream-reports/faer-029-evd-input-not-scaled`](upstream-reports/faer-029-evd-input-not-scaled/). |
 | dqds for values-only singular values | LAPACK's `dgesdd(jobz='N')` never applies a rotation: `DBDSDC(compq='N')` -> `DLASDQ` -> `DBDSQR`'s `ROTATE`-false branch -> `DLASQ1`, the dqds iteration. The fork had no such path. `bidiag_svd::dqds` now applies the shifted qd transform to the squared factors (one division and one multiplication per entry per sweep) with scaling before squaring, a positivity-checked shift retry, splitting at zero off-diagonals and Li's deflation; a rank-deficient bidiagonal keeps the rotation path. Verified against `thin_svd` for n = 2..40 over four spectra at 1e-10. It is serial, and the divide-and-conquer it replaces threads its merge, so the route is by (budget, size): dqds at one thread or below n = 256, divide and conquer above. Campaign `svd_values` ratio (below 1 is faster), 192 cells, **110 faster and none slower beyond 5%**: 32 0.909-0.911, 64 0.873-0.879, 96 0.877, 128 0.789 (1T) to 0.911, 256 0.856 (1T) and 0.999-1.004 threaded, 512 0.922 (1T) and 0.999-1.000 threaded, 1024x64 0.975-0.986, 4096x256 0.978-0.999. The 36 `qr_factor` cells the arm flagged are the binary-layout effect: a control with the kernel compiled but not called moves them by the same amount (`on/off` 1.01-1.20 while `off/pin` is 0.87-1.06). | EXPERIMENT | tlinalg-rs#46 | The predeclared target (`svd_values/64` approaching 0.275 ms) is **not** met (0.47). Three tuning attempts (shift-relative deflation, late-failure shift correction, analytic two-value deflation) each fail to converge on a wide-spectrum block: the missing piece is an admissible shift (`dlasq4`), which is also what 512+ needs. Worklog: `docs/worklogs/2026-10-23-dqds-values-only.md`; evidence `20261023T000823Z-dqds-blast` / `20261023T012045Z-dqds-blast-baseline` / `20261023T043848Z-dqds-off-control`. Tests: `dqds::tests::{matches_the_reference_on_bidiagonals,diagonal_is_its_own_answer,scale_does_not_change_the_answer}`. |
 | blocked bidiagonalisation, the dqds shift, the work copy's leading dimension, and the caller-held worker group (values-only SVD) | `bidiag_in_place` is LAPACK's `dgebrd`/`dlabrd` form: panels of 16 reduced against deferred corrections `X`/`Y`, one trailing update per panel through two plain `matmul` calls, the unblocked loop on the tail below a remaining dimension of 128 (sizes at or below 128 untouched); each panel row is gathered into a contiguous buffer, because a row of `A` has stride `m` and neither the vector products nor the reflector norm take a strided operand. dqds: ping-pong buffers instead of save/restore, release relative to the eigenvalue (`(100 eps)^2 (sigma + a)`), trailing 2x2 in closed form, the shift backed off from the trailing 2x2's eigenvalue by twice its eigenvector-weighted second-order coupling, a rejection at the last factor corrected to `tau + delta`: 3.5 accepted and 0.2 rejected sweeps per value at 512 (from 5.7 and 3.8); values-only takes dqds at every size and budget. The SVD's work copy gets 64 more bytes per column when its stride would be a multiple of 1 KiB (a 4 KiB stride put a row's entries in one L1 set: bidiag 512 8T 18.6 -> 15.3 ms). The worker group is the caller's: the library takes no `spindle::with_lock`, re-exports `faer::spindle` and documents that a caller with a thread budget holds it (the two raw-rayon loops of the block-apply split now go through spindle, as they starved under a held group). Final campaign against the pin under a harness that holds the group (`20261023T125347Z-blocked-dqds-caller-lock-2` vs `20261023T112343Z-caller-lock-baseline`, 384 cells, control 1.0027, **266 faster / 0 slower beyond 5%**): `svd_values` 32-96 0.60-0.63, 128 0.57-0.67, 256 0.60-0.67, 512 0.67 (1T) / 0.61 (8T), 4096x256 0.92 (1T) / 0.57 (8T); `svd_thin` 0.79 at 512 8T, 0.51 at 4096x256 8T. Against MKL in the same run (F order): f64 512 0.85 / 0.91 / 0.97 at 1T / 4T / 8T, 4096x256 0.83 / 0.93 / 1.11, 64 1.17; complex faster at 1T everywhere (0.60-0.99) and at 8T for 512 and both tall shapes. The unblocked loop keeps its row in a contiguous buffer too (`41161e9`: `svd_values` 64 0.296 -> 0.239 ms, 128 1.30 -> 1.01), and a basis of up to eight columns forms its `T` by inner products on the calling thread (`2351819`: the recursive QR's deepest upgrades, 2.0 -> 0.27 ms at 4096x256 8T, `qr_factor` 0.47-0.96 of the pin). Final run `20261023T193712Z-row-buffer-narrow-basis` against same-harness baselines: SVD families 363 faster / 0 slower over 384 cells (control 1.0010), `qr_factor` 174 / 0 over 192 (control 0.9947); `svd_values` 0.51 of the pin at 64, 0.61 at 512 8T, 0.56 at 4096x256 8T; against MKL f64 **64 0.94**, 512 0.85-0.95, 4096x256 0.81 / 0.92 / 1.09. Provenance: published algorithms (Dongarra-Hammarling-Sorensen 1989; Fernando-Parlett 1994; Parlett-Marques LAWN 155) and measured breakdowns only; no vendor or LAPACK code consulted, the vendor timed through its public API. | EXPERIMENT | tlinalg-rs#46 | The predeclared target (`svd_values/64` approaching 0.275 ms) is **met: 0.239 ms against the vendor's 0.255** (`41161e9`). Left: 32 (1.09), the eight-thread tall shape (1.09) and 256 at 8T (1.02), in the QR panel's narrow products and the 256-sized bidiagonalisation that does not thread (`docs/worklogs/2026-10-23-svd-64-and-tall.md`). For the maintainer: adoption of the whole set, and whether `faer::spindle` stays a public re-export. Worklogs: `docs/worklogs/2026-10-23-blocked-bidiag.md`, `docs/worklogs/2026-10-23-dqds-shift.md`; evidence in `faer-benchmark`: `20261023T091426Z-blocked-bidiag`, `20261023T100330Z-blocked-dqds` (old harness, no held group), `20261023T112343Z-caller-lock-baseline`, `20261023T120538Z-blocked-dqds-caller-lock`, `20261023T125347Z-blocked-dqds-caller-lock-2`. Tests: bidiag reconstruction at 200/300/400x300 real and complex, sequential and threaded; `dqds::tests::matches_the_reference_on_a_dense_bidiagonal` (1e-12 per value up to 200); lib 218 passed / 2 ignored. |
 | dqds falls back to the rotation path when it does not finish (#110) | Since the values-only route went to dqds (row above), a bidiagonal whose singular values span more than the range of their squares returned `NoConvergence`. The case is LAPACK's DCHKBD type 16, entries `eps^U(-2, 2)`: the zero-shift sweep's running pivot underflows to zero and every retry is rejected. It hit all four dtypes at n = 127 to 129 in `dchkbd_15_to_19` and `bidiagonal_family`. `compute_bidiag_real_svd` now finishes on the rotation path when `dqds::singular_values` fails. That is LAPACK's structure: `DBDSQR` falls through to its implicit QR when `DLASQ1` reports that dqds did not finish (`INFO = 2`). `dqds::singular_values` writes `d` only on success and never writes `e`, so the fallback starts from the caller's bidiagonal. `DLASQ1`'s scaling (the largest entry to `sqrt(eps / safmin)` before squaring) was tried and does not rescue this matrix, so it is not adopted. Fork-only: the values-only dqds route came from #52 and upstream faer has no dqds, so upstream is not affected. | KEEP | tensor4all/faer#110 | Correctness, no timing claim: only the failing inputs take the second path. Test `dchkbd_type_16_values_only_finishes` (n = 128 at recursion threshold 8 and n = 127 at 128, so both fallback routes; all four `U`/`V` modes, every dtype) runs with the regular tests and fails without the fallback. The opt-in `dchkbd_15_to_19` and `bidiagonal_family` pass in all four dtypes. Record: [worklog](docs/worklogs/2026-10-29-dqds-fallback.md). |

docs/lapack-validation.md

@@ -130,11 +130,18 @@ fixtures and with the same residual formula:
   grows like an accumulated rounding error and still catches a residual that grows with n:
   the clamp defect fixed in #121 gave 62.5 at n = 63, which a linear `20 * n / 16` (78.75)
   would have let pass. With it the f32, f64 and c64 groups and every parallel group pass.
-  What remains is c32 only: kinds 6 and 7 (scaled to `sqrt(MAX) / 16` and
-  `sqrt(MIN) * 16`) at n = 63 to 65, ratios 41.7 to 59.4 against 39.7 to 40.3.
-  Single-precision LAPACK `cgeev` (through scipy, OpenBLAS) gives 29 to 40 on the same
-  matrices, itself 39.77 against 39.69 on kind 6 at n = 63 (left vectors). That test stays
-  in `scripts/lapack-known-failures.json` under #112.
+  What remained was c32 only: kinds 6 and 7 (scaled to `sqrt(MAX) / 16` and
+  `sqrt(MIN) * 16`) at n = 63 to 65, ratios 41.7 to 59.4 against 39.7 to 40.3, where
+  single-precision LAPACK `cgeev` (through scipy, OpenBLAS) gives 29 to 40. A stage-by-stage
+  comparison with `cgeev` placed the difference in the complex Schur iteration (backward
+  error 226 to 258 against 67 to 93 on the same matrix scaled to unit size): `?GEEV` scales
+  its input into `[smlnum, bignum]` (`smlnum = sqrt(safmin) / eps`) before the Hessenberg
+  reduction, and faer did not. `evd_imp` now does the same, by a power of two, and scales the
+  eigenvalues back; the c32 group passes and the known-failures entry is gone. A sweep over
+  f32/f64/c32/c64, n = 16 to 257 and scales up to the range limits showed that the missing
+  scaling was far more serious below: on random, Grcar and near-Jordan matrices the ratio grew
+  from about `2^-90` (32-bit), and at `2^-100` (32-bit) and `2^-960` (64-bit) and below it was
+  `1e4` to `2e15` where LAPACK gives 2 to 70. Regression test: `geev_scaling_tests`.
 
 ### Phase 4 — Balancing and back-transformation