From 3c72e4327b1c892032bf57f0c0be21a514131d7d Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Thu, 29 Oct 2026 19:39:57 +0900 Subject: [PATCH 1/6] evd: scale the input into range as ?GEEV does (#112) evd_imp reduced and iterated on the matrix as given. ?GEEV first scales it into [smlnum, bignum] (smlnum = sqrt(safmin) / eps). Without that, entries below about 2^-90 (32-bit) or 2^-960 (64-bit) went through the Hessenberg reduction and the Schur iteration in subnormals and the eigenvectors lost every digit (GET22 1e5 to 2e15, LAPACK 2 to 70), and the complex Schur iteration near the top of the range lost the factor that kept the c32 LAPACK group failing. The scaling is a power of two, so it is exact, and the eigenvalues are scaled back. Co-Authored-By: Claude Opus 5.5 --- docs/lapack-validation.md | 17 ++-- faer/src/linalg/evd/mod.rs | 140 +++++++++++++++++++++++++++++ scripts/lapack-known-failures.json | 10 +-- 3 files changed, 153 insertions(+), 14 deletions(-) diff --git a/docs/lapack-validation.md b/docs/lapack-validation.md index caa6e546..a683b9e8 100644 --- a/docs/lapack-validation.md +++ b/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: every eigenvector digit was lost for entries under + about `2^-90` (32-bit) and `2^-960` (64-bit), GET22 ratios of `1e5` to `2e15` where LAPACK + gives 2 to 70. Regression test: `geev_scaling_tests`. ### Phase 4 — Balancing and back-transformation diff --git a/faer/src/linalg/evd/mod.rs b/faer/src/linalg/evd/mod.rs index 7aba8438..6f51ad8a 100644 --- a/faer/src/linalg/evd/mod.rs +++ b/faer/src/linalg/evd/mod.rs @@ -1833,6 +1833,44 @@ pub fn eigenvectors_from_schur_cplx( Err(EvdError::NoConvergence) } } +/// `?GEEV`'s scaling of its input: when the largest entry modulus is below +/// `smlnum = sqrt(min_positive) / eps` or above `bignum = 1 / smlnum`, the power of two that +/// brings it to just inside that bound, else `None`. LAPACK scales by `smlnum / anrm` (or +/// `bignum / anrm`); a power of two lands within a factor of two of the same bound and scales +/// every entry exactly. Without it the Hessenberg reduction and the Schur iteration of a matrix +/// whose entries are near the bottom of the range work in subnormals, and the eigenvectors lost +/// every digit (a GET22 ratio of about `1 / eps` below `2^-90` for `f32` and `2^-960` for +/// `f64`); near the top, the complex Schur iteration lost a factor of three against `cgeev` +/// (#112). +fn geev_scaling(A: MatRef<'_, T>) -> Option { + let n = A.nrows(); + let mut anrm = zero::(); + for j in 0..n { + for i in 0..n { + anrm = anrm.fmax(A[(i, j)].abs()); + } + } + let smlnum = min_positive::().sqrt() / eps::(); + let bignum = smlnum.recip(); + let two = from_f64::(2.0); + let half = from_f64::(0.5); + if anrm == zero::() || (anrm >= smlnum && anrm <= bignum) { + return None; + } + let mut f = one::(); + if anrm < smlnum { + // `anrm * f` in `[smlnum, 2 smlnum)` + while &anrm * &f < smlnum { + f = f * &two; + } + } else { + // `anrm * f` in `(bignum / 2, bignum]` + while &anrm * &f > bignum { + f = f * ½ + } + } + Some(f) +} fn evd_imp( A: MatRef<'_, T>, s: ColMut<'_, T>, @@ -1879,6 +1917,11 @@ fn evd_imp( None }; H.copy_from(A); + // the eigenvectors of `f A` are those of `A`; its eigenvalues are scaled back at the end + let scaling = geev_scaling(A); + if let Some(f) = &scaling { + z!(H.rb_mut()).for_each(|uz!(h)| *h = h.mul_real(f)); + } { let (mut householder, stack) = unsafe { temp_mat_uninit::(bs, n - 1, stack) }; @@ -2013,6 +2056,13 @@ fn evd_imp( return Err(EvdError::NoConvergence); } } + if let Some(f) = &scaling { + let f = f.recip(); + z!(s.rb_mut()).for_each(|uz!(x)| *x = x.mul_real(&f)); + if let Some(s_im) = s_im.rb_mut() { + z!(s_im).for_each(|uz!(x)| *x = x.mul_real(&f)); + } + } Ok(()) } /// computes the matrix $A$'s eigendecomposition @@ -3101,6 +3151,96 @@ mod tridiagonal_range_tests { } } +#[cfg(test)] +mod geev_scaling_tests { + //! Regression tests for the input scaling of `evd_imp` (#112). A matrix whose largest entry is + //! near the bottom of the range went through the Hessenberg reduction and the Schur iteration + //! in subnormals: a Grcar or random matrix scaled by `2^-120` (`f32`, `c32`) or `2^-1010` + //! (`f64`, `c64`) had a GET22 ratio of about `1 / eps` (`4e6` and `2e15`). Near the top the + //! complex Schur iteration lost a factor of three against `cgeev`. Residuals are measured in + //! `f64` on the dimensionless matrix `a / scale`, exact for a power of two. + use crate::{assert, c32, c64, Mat, Par}; + use num_complex::Complex64 as Z; + + /// GET22's `||a u - u diag(w)||_1 / (||a||_1 ||u||_1 eps)` for one side + fn get22(a: &[Vec], u: &[Vec], w: &[Z], eps: f64) -> f64 { + let n = a.len(); + let (mut an, mut un, mut err) = (0.0f64, 0.0f64, 0.0f64); + for j in 0..n { + let (mut ac, mut uc, mut ec) = (0.0, 0.0, 0.0); + for i in 0..n { + ac += a[i][j].norm(); + uc += u[i][j].norm(); + let mut r = -u[i][j] * w[j]; + for k in 0..n { + r += a[i][k] * u[k][j]; + } + ec += r.norm(); + } + an = an.max(ac); + un = un.max(uc); + err = err.max(ec); + } + err / an / un / eps + } + + /// the Grcar matrix (ones on the diagonal and the three superdiagonals, minus one on the + /// subdiagonal) and a pseudo-random matrix, both with largest entry one + fn unit(n: usize, cplx: bool) -> [Vec>; 2] { + let grcar = (0..n) + .map(|i| { + (0..n) + .map(|j| match j as isize - i as isize { + -1 => Z::new(-1.0, 0.0), + 0..=3 => Z::new(1.0, 0.0), + _ => Z::new(0.0, 0.0), + }) + .collect() + }) + .collect(); + let mut state = 0x2545_f491_4f6c_dd1du64; + let mut next = move || { + state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + // a multiple of `2^-10` in `[-1, 1)`, exact in every precision + ((state >> 54) as f64 - 512.0) / 512.0 + }; + let rand = (0..n) + .map(|_| (0..n).map(|_| Z::new(next(), if cplx { next() } else { 0.0 })).collect()) + .collect(); + [grcar, rand] + } + + macro_rules! check { + ($t:ty, $mk:expr, $z:expr, $eps:expr, $scales:expr, $cplx:expr) => {{ + let n = 32; + // the GET22 threshold of the LAPACK nonsymmetric EVD tests + let thresh = 20.0 * f64::max(1.0, (n as f64 / 16.0).sqrt()); + for unit in unit(n, $cplx) { + for e in $scales { + let scale = 2.0f64.powi(e); + let a = Mat::<$t>::from_fn(n, n, |i, j| $mk(unit[i][j] * scale)); + for (side, m) in [("right", a.clone()), ("left", a.adjoint().to_owned())] { + let dimless: Vec> = (0..n).map(|i| (0..n).map(|j| $z(m[(i, j)]) / scale).collect()).collect(); + let evd = m.eigen(Par::Seq).unwrap(); + let w: Vec = (0..n).map(|i| Z::new(evd.S()[i].re as f64, evd.S()[i].im as f64) / scale).collect(); + let u: Vec> = (0..n).map(|i| (0..n).map(|j| Z::new(evd.U()[(i, j)].re as f64, evd.U()[(i, j)].im as f64)).collect()).collect(); + let r = get22(&dimless, &u, &w, $eps); + assert!(r <= thresh, "{} scale 2^{e} {side}: GET22 {r}", stringify!($t)); + } + } + } + }}; + } + + #[test] + fn evd_scales_its_input_into_range() { + check!(f32, |z: Z| z.re as f32, |x: f32| Z::new(x as f64, 0.0), f32::EPSILON as f64, [-120, -100, 0, 100, 120], false); + check!(c32, |z: Z| c32::new(z.re as f32, z.im as f32), |x: c32| Z::new(x.re as f64, x.im as f64), f32::EPSILON as f64, [-120, -100, 0, 100, 120], true); + check!(f64, |z: Z| z.re, |x: f64| Z::new(x, 0.0), f64::EPSILON, [-1010, -1000, -960, 0, 960, 1015], false); + check!(c64, |z: Z| c64::new(z.re, z.im), |x: c64| Z::new(x.re, x.im), f64::EPSILON, [-1010, -1000, -960, 0, 960, 1015], true); + } +} + #[cfg(test)] #[path = "tests/mod.rs"] mod lapack; diff --git a/scripts/lapack-known-failures.json b/scripts/lapack-known-failures.json index c0bafe90..d3cff47b 100644 --- a/scripts/lapack-known-failures.json +++ b/scripts/lapack-known-failures.json @@ -5,13 +5,5 @@ "failure that is not listed here, or for an entry that has started passing (so a fix has to", "remove its entry in the same change). Never add a test without the issue that describes it." ], - "known": [ - { - "issue": 112, - "note": "c32 only, after the size-dependent GET22 threshold: kinds 6 and 7 (scaled to sqrt(MAX)/16 and sqrt(MIN)*16) at n = 63-65 give ratios of 41.7-59.4 against 39.7-40.3; single-precision LAPACK cgeev gives 29-40 on the same matrices", - "tests": [ - "linalg::evd::lapack::c32::lapack_nonsymmetric_evd_cases" - ] - } - ] + "known": [] } From 798e81e6ca7cb173856516bb77be870b89a3a43e Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Thu, 29 Oct 2026 19:40:57 +0900 Subject: [PATCH 2/6] upstream-reports, PATCHSET: the EVD input scaling (#112) faer-029 reproduces the lost eigenvectors of a tiny matrix on faer 0.24.4 (f64 2^-960: GET22 1.1e15; f32 2^-100: 1.8e6; LAPACK 16 and 16) and on the fork prints 17.7 and 15.2 at every scale. Co-Authored-By: Claude Opus 5.5 --- PATCHSET.md | 1 + upstream-reports/README.md | 1 + .../faer-029-evd-input-not-scaled/Cargo.lock | 678 ++++++++++++++++++ .../faer-029-evd-input-not-scaled/Cargo.toml | 11 + .../faer-029-evd-input-not-scaled/README.md | 52 ++ .../faer-029-evd-input-not-scaled/src/main.rs | 50 ++ 6 files changed, 793 insertions(+) create mode 100644 upstream-reports/faer-029-evd-input-not-scaled/Cargo.lock create mode 100644 upstream-reports/faer-029-evd-input-not-scaled/Cargo.toml create mode 100644 upstream-reports/faer-029-evd-input-not-scaled/README.md create mode 100644 upstream-reports/faer-029-evd-input-not-scaled/src/main.rs diff --git a/PATCHSET.md b/PATCHSET.md index 49c85612..65b0bea5 100644 --- a/PATCHSET.md +++ b/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::()` 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 below about `2^-90` (32-bit) or `2^-960` (64-bit) was reduced and iterated in subnormals and every eigenvector digit was lost (GET22 `1e5` to `2e15`, LAPACK 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 (exact) and scales the eigenvalues back. | 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 and left vectors, GET22 below `20 * max(1, sqrt(n / 16))`) fails without the scaling (`f32`, `2^-120`: 1.7e6) and passes with it. 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). | diff --git a/upstream-reports/README.md b/upstream-reports/README.md index 75eed75a..7224376c 100644 --- a/upstream-reports/README.md +++ b/upstream-reports/README.md @@ -42,6 +42,7 @@ The program prints the expected and the actual value, and exits with status 1 wh | [faer-026](faer-026-svd-qr-two-by-two-no-convergence/) | The SVD returns `NoConvergence` on a 12 x 10 complex matrix | `faer` 0.24.4 | no convergence | tensor4all/faer#84 (`f83ffb0`) | not yet | 2026-10-29 | | [faer-027](faer-027-qr-near-dependent-column-skip/) | Unpivoted QR skips a nearly dependent column of a tall matrix | `faer` 0.24.4 | inaccurate value | tensor4all/faer#1 (`2e6ea52`, `d16929b`) | not yet | 2026-10-29 | | [faer-028](faer-028-sparse-lu-numeric-singular/) | Sparse LU factors a numerically singular matrix without an error | `faer` 0.24.4 | NaN | tensor4all/faer#132 | not yet | 2026-10-29 | +| [faer-029](faer-029-evd-input-not-scaled/) | The nonsymmetric EVD loses every eigenvector digit when the entries are tiny | `faer` 0.24.4 | inaccurate value | tensor4all/faer#112 | not yet | 2026-10-29 | | [gemm-x86-001](gemm-x86-001-diag-negative-stride/) | A product with a diagonal gives a wrong value when the left matrix is stored backwards | `private-gemm-x86` 0.1.22 | wrong value | tensor4all/t4a-private-gemm-x86#5 (`84ca6db`) | not yet | 2026-10-29 | "Fixed upstream" was checked against upstream `main`: `faer` at `0ecd846` (2026-10-22) and `private-gemm-x86` at `71872b9` (its 0.1.22 release). diff --git a/upstream-reports/faer-029-evd-input-not-scaled/Cargo.lock b/upstream-reports/faer-029-evd-input-not-scaled/Cargo.lock new file mode 100644 index 00000000..a0a7f71e --- /dev/null +++ b/upstream-reports/faer-029-evd-input-not-scaled/Cargo.lock @@ -0,0 +1,678 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 4 + +[[package]] +name = "autocfg" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" + +[[package]] +name = "bitflags" +version = "2.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3ded4057c258ba199e2d26386d3af3780957ecaee6c4ef4041c6b4b8b97c0b06" + +[[package]] +name = "bytemuck" +version = "1.25.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" +dependencies = [ + "bytemuck_derive", +] + +[[package]] +name = "bytemuck_derive" +version = "1.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6a1f896587b6f2c069c73d2f0913e2d590c3990285cd2f0b6aa02b786b4c679c" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.9", +] + +[[package]] +name = "byteorder" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" + +[[package]] +name = "cfg-if" +version = "1.0.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4e7648175b45a9a48536d676f68d918270699102aa8dab5496df06904c914600" + +[[package]] +name = "crunchy" +version = "0.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "460fbee9c2c2f33933d720630a6a0bac33ba7053db5344fac858d4b8952d77d5" + +[[package]] +name = "defer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "930c7171c8df9fb1782bdf9b918ed9ed2d33d1d22300abb754f9085bc48bf8e8" + +[[package]] +name = "dyn-stack" +version = "0.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1c4713e43e2886ba72b8271aa66c93d722116acf7a75555cce11dcde84388fe8" +dependencies = [ + "bytemuck", + "dyn-stack-macros", +] + +[[package]] +name = "dyn-stack-macros" +version = "0.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e1d926b4d407d372f141f93bb444696142c29d32962ccbd3531117cf3aa0bfa9" + +[[package]] +name = "enum-as-inner" +version = "0.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a1e6a265c649f3f5979b601d26f1d05ada116434c87741c9493cb56218f76cbc" +dependencies = [ + "heck", + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "equator" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c35da53b5a021d2484a7cc49b2ac7f2d840f8236a286f84202369bd338d761ea" +dependencies = [ + "equator-macro 0.2.1", +] + +[[package]] +name = "equator" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "02da895aab06bbebefb6b2595f6d637b18c9ff629b4cd840965bb3164e4194b0" +dependencies = [ + "equator-macro 0.6.0", +] + +[[package]] +name = "equator-macro" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3bf679796c0322556351f287a51b49e48f7c4986e727b5dd78c972d30e2e16cc" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "equator-macro" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b14b339eb76d07f052cdbad76ca7c1310e56173a138095d3bf42a23c06ef5d8" + +[[package]] +name = "faer" +version = "0.24.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5ab6df3dd147fe8d702a288b95bcd8fcc499ab572fc80da6828f60cd4d524d67" +dependencies = [ + "bytemuck", + "dyn-stack", + "equator 0.6.0", + "faer-traits", + "gemm", + "generativity", + "libm", + "nano-gemm", + "num-complex", + "num-traits", + "private-gemm-x86", + "pulp", + "reborrow", +] + +[[package]] +name = "faer-029-evd-input-not-scaled" +version = "0.0.0" +dependencies = [ + "faer", +] + +[[package]] +name = "faer-traits" +version = "0.24.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b87d23ed7ab1f26c0cba0e5b9e061a796fbb7dc170fa8bee6970055a1308bb0f" +dependencies = [ + "bytemuck", + "dyn-stack", + "generativity", + "libm", + "num-complex", + "num-traits", + "pulp", + "qd", + "reborrow", +] + +[[package]] +name = "gemm" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa0673db364b12263d103b68337a68fbecc541d6f6b61ba72fe438654709eacb" +dependencies = [ + "dyn-stack", + "gemm-c32", + "gemm-c64", + "gemm-common", + "gemm-f16", + "gemm-f32", + "gemm-f64", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-c32" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "086936dbdcb99e37aad81d320f98f670e53c1e55a98bee70573e83f95beb128c" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-c64" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "20c8aeeeec425959bda4d9827664029ba1501a90a0d1e6228e48bef741db3a3f" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-common" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "88027625910cc9b1085aaaa1c4bc46bb3a36aad323452b33c25b5e4e7c8e2a3e" +dependencies = [ + "bytemuck", + "dyn-stack", + "half", + "libm", + "num-complex", + "num-traits", + "once_cell", + "paste", + "pulp", + "raw-cpuid", + "seq-macro", + "sysctl", +] + +[[package]] +name = "gemm-f16" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e3df7a55202e6cd6739d82ae3399c8e0c7e1402859b30e4cb780e61525d9486e" +dependencies = [ + "dyn-stack", + "gemm-common", + "gemm-f32", + "half", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-f32" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "02e0b8c9da1fbec6e3e3ab2ce6bc259ef18eb5f6f0d3e4edf54b75f9fd41a81c" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-f64" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "056131e8f2a521bfab322f804ccd652520c79700d81209e9d9275bbdecaadc6a" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "generativity" +version = "1.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d2c81fb5260e37854d09d5c87183309fd8c555b75289427884b25660bc87a85e" + +[[package]] +name = "half" +version = "2.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6ea2d84b969582b4b1864a92dc5d27cd2b77b622a8d79306834f1be5ba20d84b" +dependencies = [ + "bytemuck", + "cfg-if", + "crunchy", + "num-traits", + "zerocopy", +] + +[[package]] +name = "heck" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" + +[[package]] +name = "hermit-abi" +version = "0.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e17592d60ebacc7d5e169f4663c5f84f9161cc90328abcfe8456f41e4dfcb284" + +[[package]] +name = "interpol" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eb58032ba748f4010d15912a1855a8a0b1ba9eaad3395b0c171c09b3b356ae50" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "libc" +version = "0.2.191" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "65cd87c9a6d60536580b316ad985726104730a5151c99845f2490428ef91240a" + +[[package]] +name = "libm" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6d2cec3eae94f9f509c767b45932f1ada8350c4bdb85af2fcab4a3c14807981" + +[[package]] +name = "nano-gemm" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9e04345dc84b498ff89fe0d38543d1f170da9e43a2c2bcee73a0f9069f72d081" +dependencies = [ + "equator 0.2.2", + "nano-gemm-c32", + "nano-gemm-c64", + "nano-gemm-codegen", + "nano-gemm-core", + "nano-gemm-f32", + "nano-gemm-f64", + "num-complex", +] + +[[package]] +name = "nano-gemm-c32" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0775b1e2520e64deee8fc78b7732e3091fb7585017c0b0f9f4b451757bbbc562" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", + "num-complex", +] + +[[package]] +name = "nano-gemm-c64" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9af49a20d58816e6b5ee65f64142e50edb5eba152678d4bb7377fcbf63f8437a" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", + "num-complex", +] + +[[package]] +name = "nano-gemm-codegen" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6cc8d495c791627779477a2cf5df60049f5b165342610eb0d76bee5ff5c5d74c" + +[[package]] +name = "nano-gemm-core" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d998dfa644de87a0f8660e5ea511d7cb5c33b5a2d9847b7af57a2565105089f0" + +[[package]] +name = "nano-gemm-f32" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "879d962e79bc8952e4ad21ca4845a21132540ed3f5e01184b2ff7f720e666523" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", +] + +[[package]] +name = "nano-gemm-f64" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9a513473dce7dc00c7e7c318481ca4494034e76997218d8dad51bd9f007a815" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", +] + +[[package]] +name = "num-complex" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" +dependencies = [ + "bytemuck", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", + "libm", +] + +[[package]] +name = "num_cpus" +version = "1.17.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91df4bbde75afed763b708b7eee1e8e7651e02d97f6d5dd763e89367e957b23b" +dependencies = [ + "hermit-abi", + "libc", +] + +[[package]] +name = "once_cell" +version = "1.21.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50" + +[[package]] +name = "paste" +version = "1.0.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" + +[[package]] +name = "private-gemm-x86" +version = "0.1.22" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bba0a7831d1297712e7e7b33759b0c837941e79196dd05df506d5f36f7f78619" +dependencies = [ + "defer", + "interpol", + "num_cpus", + "raw-cpuid", +] + +[[package]] +name = "proc-macro2" +version = "1.0.107" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "985e7ec9bb745e6ce6535b544d84d6cd6f7ad8bd711c398938ae983b91a766d9" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "pulp" +version = "0.22.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "046aa45b989642ec2e4717c8e72d677b13edd831a4d3b6cf37d9a3e54912496a" +dependencies = [ + "bytemuck", + "cfg-if", + "libm", + "num-complex", + "paste", + "pulp-wasm-simd-flag", + "raw-cpuid", + "reborrow", + "version_check", +] + +[[package]] +name = "pulp-wasm-simd-flag" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1d8f70e07b9c3962945a74e59ca1c511bba65b6419468acc217c457d93f3c740" + +[[package]] +name = "qd" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "15f1304a5aecdcfe9ee72fbba90aa37b3aa067a69d14cb7f3d9deada0be7c07c" +dependencies = [ + "bytemuck", + "libm", + "num-traits", + "pulp", +] + +[[package]] +name = "quote" +version = "1.0.47" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fbf4db142a473a8d80c26bbf18454ed458bf8d26c8219c331daecfdbd079001" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "raw-cpuid" +version = "11.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "498cd0dc59d73224351ee52a95fee0f1a617a2eae0e7d9d720cc622c73a54186" +dependencies = [ + "bitflags", +] + +[[package]] +name = "reborrow" +version = "0.5.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "03251193000f4bd3b042892be858ee50e8b3719f2b08e5833ac4353724632430" + +[[package]] +name = "same-file" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93fc1dc3aaa9bfed95e02e6eadabb4baf7e3078b0bd1b4d7b6b0b68378900502" +dependencies = [ + "winapi-util", +] + +[[package]] +name = "seq-macro" +version = "0.3.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1bc711410fbe7399f390ca1c3b60ad0f53f80e95c5eb935e52268a0e2cd49acc" + +[[package]] +name = "syn" +version = "1.0.109" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "2.0.119" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "872831b642d1a07999a962a351ed35b955ea2cfc8f3862091e2a240a84f17297" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "3.0.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d78c8dee4c7bf0e14673097256fed6142ce9d3b85a408189d07482442145823b" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "sysctl" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "01198a2debb237c62b6826ec7081082d951f46dbb64b0e8c7649a452230d1dfc" +dependencies = [ + "bitflags", + "byteorder", + "enum-as-inner", + "libc", + "thiserror", + "walkdir", +] + +[[package]] +name = "thiserror" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6aaf5339b578ea85b50e080feb250a3e8ae8cfcdff9a461c9ec2904bc923f52" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "unicode-ident" +version = "1.0.27" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a2c754d6c33795a1c324727428e5a7dedb5b06195f9890bdbcba760d3e246563" + +[[package]] +name = "version_check" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" + +[[package]] +name = "walkdir" +version = "2.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "29790946404f91d9c5d06f9874efddea1dc06c5efe94541a7d6863108e3a5e4b" +dependencies = [ + "same-file", + "winapi-util", +] + +[[package]] +name = "winapi-util" +version = "0.1.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c2a7b1c03c876122aa43f3020e6c3c3ee5c05081c9a00739faf7503aeba10d22" +dependencies = [ + "windows-sys", +] + +[[package]] +name = "windows-link" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f0805222e57f7521d6a62e36fa9163bc891acd422f971defe97d64e70d0a4fe5" + +[[package]] +name = "windows-sys" +version = "0.61.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ae137229bcbd6cdf0f7b80a31df61766145077ddf49416a728b02cb3921ff3fc" +dependencies = [ + "windows-link", +] + +[[package]] +name = "zerocopy" +version = "0.8.63" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e5fe1f8f1b06191a00962174c61aa5005e0bb391a6d80d07e24d115c01a92ed8" +dependencies = [ + "zerocopy-derive", +] + +[[package]] +name = "zerocopy-derive" +version = "0.8.63" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "863ad3ac83293fb4d740aedbfdc9240dd8d1a50c1099acd76ce80ce7c7230c7f" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] diff --git a/upstream-reports/faer-029-evd-input-not-scaled/Cargo.toml b/upstream-reports/faer-029-evd-input-not-scaled/Cargo.toml new file mode 100644 index 00000000..1b552762 --- /dev/null +++ b/upstream-reports/faer-029-evd-input-not-scaled/Cargo.toml @@ -0,0 +1,11 @@ +[package] +name = "faer-029-evd-input-not-scaled" +version = "0.0.0" +edition = "2021" +publish = false + +# Not part of the t4a-faer workspace: this reproduces the published crate on its own. +[workspace] + +[dependencies] +faer = { version = "=0.24.4", default-features = false, features = ["std", "linalg"] } diff --git a/upstream-reports/faer-029-evd-input-not-scaled/README.md b/upstream-reports/faer-029-evd-input-not-scaled/README.md new file mode 100644 index 00000000..37ae9a8e --- /dev/null +++ b/upstream-reports/faer-029-evd-input-not-scaled/README.md @@ -0,0 +1,52 @@ +# faer-029: the nonsymmetric EVD loses every eigenvector digit when the entries are tiny + +**What happens:** `Mat::eigen` returns eigenvectors whose residual `|A U - U S|` is about `1 / eps` times larger than LAPACK's, so they have no correct digit. +**Which input:** any matrix whose largest entry is below about `2^-90` (`f32`, `c32`) or `2^-960` (`f64`, `c64`). The program uses the 64 x 64 Grcar matrix scaled by a power of two. +**Which version:** `faer` 0.24.4. +**Added to this list:** 2026-10-29. + +## Reproduce + +```sh +cargo run +``` + +The program prints LAPACK's GET22 residual ratio, `|A U - U S|_1 / (|A|_1 |U|_1 eps)`, measured in `f64` on the matrix divided by its scale (exact for a power of two). +LAPACK's CHKHS test accepts a ratio below 20. + +| type | scale | expected (LAPACK `?geev`, scipy/OpenBLAS) | `faer` 0.24.4 | +|---|---|---|---| +| f64 | 1 | 16.4 | 10.6 | +| f64 | `2^-960` | 16.4 | **1.1e15** | +| f64 | `2^-1000` | 16.4 | **9.0e14** | +| f32 | 1 | 15.8 | 10.4 | +| f32 | `2^-100` | 15.8 | **1.8e6** | +| f32 | `2^-110` | 15.8 | **1.7e6** | + +A sweep over `f32`, `f64`, `c32` and `c64`, n = 16 to 257, random, Grcar, near-Jordan and graded matrices, gives the same picture: from about `2^-90` (32-bit) and `2^-960` (64-bit) down, every family loses all digits, while LAPACK stays at 2 to 70. +Near the top of the range the loss is smaller: on LAPACK's own `c32` test matrices scaled to `sqrt(MAX) / 16`, `faer` gives a ratio of 42 to 59 where `cgeev` gives 29 to 40. + +## Cause + +`evd_imp` (`src/linalg/evd/mod.rs`) copies the matrix and runs the Hessenberg reduction and the Schur iteration on it as given. +LAPACK's `?GEEV` first computes the largest entry modulus `anrm` and, when it is below `smlnum = sqrt(safmin) / eps` or above `bignum = 1 / smlnum`, scales the matrix to that bound; it scales the eigenvalues back at the end. +Without that step the reduction and the iteration of a tiny matrix work in subnormals, where the relative precision is gone. +Measured stage by stage on a `c32` test matrix near the top, the Hessenberg backward error is the same with and without scaling (16 to 20), and the Schur backward error is 226 to 258 without and 67 to 93 with it. + +## How the fork fixed it + +`evd_imp` scales its input as `?GEEV` does, with a power of two instead of `smlnum / anrm`, so the scaling is exact. +The eigenvectors do not change; the eigenvalues are multiplied back. +With that change this program prints 17.7 (`f64`) and 15.2 (`f32`) at every scale. + +Fork regression test: `geev_scaling_tests::evd_scales_its_input_into_range` in `faer/src/linalg/evd/mod.rs` (Grcar and pseudo-random matrices, all four types, scales to both ends of the range). + +## Affected versions + +Confirmed by running this program: `faer` 0.24.4. +Not fixed upstream: upstream `main` at `0ecd846` (2026-10-22) has the same `evd_imp` without input scaling. + +## References + +* Fork fix: tensor4all/faer#112. +* Found by: the bounded LAPACK validation's nonsymmetric EVD family (`c32`) and the stage-by-stage comparison with `cgeev` asked for in tensor4all/faer#112. diff --git a/upstream-reports/faer-029-evd-input-not-scaled/src/main.rs b/upstream-reports/faer-029-evd-input-not-scaled/src/main.rs new file mode 100644 index 00000000..a88caaeb --- /dev/null +++ b/upstream-reports/faer-029-evd-input-not-scaled/src/main.rs @@ -0,0 +1,50 @@ +//! faer-029: `Mat::eigen` loses every digit of the eigenvectors when the entries are tiny. +use faer::{c64, Mat}; + +/// LAPACK's GET22 residual `|A U - U S|_1 / (|A|_1 |U|_1 eps)`, in `f64`, on the dimensionless +/// matrix `a / scale` (exact for a power of two) +fn get22(a: &Mat, u: &Mat, s: &[c64], eps: f64) -> f64 { + let n = a.nrows(); + let norm1 = |m: &Mat| { + (0..n).map(|j| (0..n).map(|i| m[(i, j)].norm()).sum::()).fold(0.0, f64::max) + }; + let us = Mat::::from_fn(n, n, |i, j| u[(i, j)] * s[j]); + let residual = a * u - us; + norm1(&residual) / norm1(a) / norm1(u) / eps +} + +fn main() { + let threshold = 20.0; + let n = 64; + // the Grcar matrix: ones on the diagonal and the three superdiagonals, minus one below + let grcar = |i: usize, j: usize| match j as isize - i as isize { + -1 => -1.0, + 0..=3 => 1.0, + _ => 0.0, + }; + let mut failed = false; + println!("{:>5} | {:>6} | {:>10} | {:>12}", "type", "scale", "expected", "actual"); + for (ty, e) in [("f64", 0), ("f64", -960), ("f64", -1000), ("f32", 0), ("f32", -100), ("f32", -110)] { + let scale = 2.0f64.powi(e); + let dimless = Mat::::from_fn(n, n, |i, j| c64::new(grcar(i, j), 0.0)); + let (u, s, eps) = if ty == "f64" { + let a = Mat::::from_fn(n, n, |i, j| grcar(i, j) * scale); + let eig = a.eigen().unwrap(); + let s: Vec = (0..n).map(|i| eig.S()[i] / scale).collect(); + (eig.U().to_owned(), s, f64::EPSILON) + } else { + let a = Mat::::from_fn(n, n, |i, j| (grcar(i, j) * scale) as f32); + let eig = a.eigen().unwrap(); + let s: Vec = (0..n).map(|i| c64::new(eig.S()[i].re as f64, eig.S()[i].im as f64) / scale).collect(); + let u = Mat::::from_fn(n, n, |i, j| c64::new(eig.U()[(i, j)].re as f64, eig.U()[(i, j)].im as f64)); + (u, s, f32::EPSILON as f64) + }; + let ratio = get22(&dimless, &u, &s, eps); + println!("{ty:>5} | 2^{e:<4} | {:>10} | {ratio:>12.4e}", format!("< {threshold}")); + failed |= !(ratio < threshold); + } + if failed { + println!("the eigenvectors of a matrix with tiny entries have lost every digit"); + std::process::exit(1); + } +} From 6ac62a18daade056265651a332788064cf8da553 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Thu, 29 Oct 2026 19:54:38 +0900 Subject: [PATCH 3/6] evd: pick the input scale from the largest |re| or |im| (#112 post-review) faer's complex modulus is inf for 2^1023 + 0i and zero for 2^-1074 + 0i, so a modulus-based scan selected a zero factor (NaN eigenvalue with Ok) or no scaling at all. The tests reject nonfinite results before f64::max can drop them; the report notes that faer 0.24.4 returns nonfinite eigenvectors with Ok at 2^-1000 (f64) and 2^-100 (f32); the docs restrict their claims to the measured families. Co-Authored-By: Claude Opus 5.5 --- PATCHSET.md | 2 +- docs/lapack-validation.md | 6 +- faer/src/linalg/evd/mod.rs | 59 ++++++++++++++++--- .../faer-029-evd-input-not-scaled/README.md | 18 +++--- .../faer-029-evd-input-not-scaled/src/main.rs | 4 +- 5 files changed, 66 insertions(+), 23 deletions(-) diff --git a/PATCHSET.md b/PATCHSET.md index 65b0bea5..6bf94f7d 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -101,7 +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::()` 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 below about `2^-90` (32-bit) or `2^-960` (64-bit) was reduced and iterated in subnormals and every eigenvector digit was lost (GET22 `1e5` to `2e15`, LAPACK 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 (exact) and scales the eigenvalues back. | 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 and left vectors, GET22 below `20 * max(1, sqrt(n / 16))`) fails without the scaling (`f32`, `2^-120`: 1.7e6) and passes with it. 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/). | +| 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). | diff --git a/docs/lapack-validation.md b/docs/lapack-validation.md index a683b9e8..0a235151 100644 --- a/docs/lapack-validation.md +++ b/docs/lapack-validation.md @@ -139,9 +139,9 @@ fixtures and with the same residual formula: 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: every eigenvector digit was lost for entries under - about `2^-90` (32-bit) and `2^-960` (64-bit), GET22 ratios of `1e5` to `2e15` where LAPACK - gives 2 to 70. Regression test: `geev_scaling_tests`. + 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 diff --git a/faer/src/linalg/evd/mod.rs b/faer/src/linalg/evd/mod.rs index 6f51ad8a..8347365f 100644 --- a/faer/src/linalg/evd/mod.rs +++ b/faer/src/linalg/evd/mod.rs @@ -1833,21 +1833,24 @@ pub fn eigenvectors_from_schur_cplx( Err(EvdError::NoConvergence) } } -/// `?GEEV`'s scaling of its input: when the largest entry modulus is below +/// `?GEEV`'s scaling of its input: when the largest entry is below /// `smlnum = sqrt(min_positive) / eps` or above `bignum = 1 / smlnum`, the power of two that /// brings it to just inside that bound, else `None`. LAPACK scales by `smlnum / anrm` (or -/// `bignum / anrm`); a power of two lands within a factor of two of the same bound and scales -/// every entry exactly. Without it the Hessenberg reduction and the Schur iteration of a matrix -/// whose entries are near the bottom of the range work in subnormals, and the eigenvectors lost -/// every digit (a GET22 ratio of about `1 / eps` below `2^-90` for `f32` and `2^-960` for -/// `f64`); near the top, the complex Schur iteration lost a factor of three against `cgeev` -/// (#112). +/// `bignum / anrm`); a power of two lands within a factor of two of the same bound and does not +/// round (an entry far below the largest can still underflow when scaling down, as in LAPACK). +/// The largest entry is the largest `|re|` or `|im|`, within a factor `sqrt(2)` of LAPACK's +/// modulus: faer's complex modulus is itself inf for `2^1023 + 0i` and zero for `2^-1074 + 0i`, +/// which would select a zero factor or none. Without the scaling, the Hessenberg reduction and +/// the Schur iteration of a matrix whose entries are near the bottom of the range work in +/// subnormals: on Grcar, random and near-Jordan matrices the GET22 ratio grew from about `2^-90` +/// (32-bit) and reached `1 / eps` by `2^-100` (32-bit) and `2^-960` (64-bit); near the top, the +/// complex Schur iteration lost a factor of three against `cgeev` (#112). fn geev_scaling(A: MatRef<'_, T>) -> Option { let n = A.nrows(); let mut anrm = zero::(); for j in 0..n { for i in 0..n { - anrm = anrm.fmax(A[(i, j)].abs()); + anrm = anrm.fmax(A[(i, j)].real().abs()).fmax(A[(i, j)].imag().abs()); } } let smlnum = min_positive::().sqrt() / eps::(); @@ -3158,7 +3161,9 @@ mod geev_scaling_tests { //! in subnormals: a Grcar or random matrix scaled by `2^-120` (`f32`, `c32`) or `2^-1010` //! (`f64`, `c64`) had a GET22 ratio of about `1 / eps` (`4e6` and `2e15`). Near the top the //! complex Schur iteration lost a factor of three against `cgeev`. Residuals are measured in - //! `f64` on the dimensionless matrix `a / scale`, exact for a power of two. + //! `f64` on the dimensionless matrix `a / scale`, exact for a power of two. The left side is + //! checked as the right eigenvectors of `a^H`; the `u_left` branch of `evd_imp` itself is + //! exercised by the LAPACK nonsymmetric EVD tests. use crate::{assert, c32, c64, Mat, Par}; use num_complex::Complex64 as Z; @@ -3224,6 +3229,8 @@ mod geev_scaling_tests { let evd = m.eigen(Par::Seq).unwrap(); let w: Vec = (0..n).map(|i| Z::new(evd.S()[i].re as f64, evd.S()[i].im as f64) / scale).collect(); let u: Vec> = (0..n).map(|i| (0..n).map(|j| Z::new(evd.U()[(i, j)].re as f64, evd.U()[(i, j)].im as f64)).collect()).collect(); + // `f64::max` drops a NaN, so a nonfinite result is rejected first + assert!(w.iter().all(|w| w.is_finite()) && u.iter().flatten().all(|u| u.is_finite())); let r = get22(&dimless, &u, &w, $eps); assert!(r <= thresh, "{} scale 2^{e} {side}: GET22 {r}", stringify!($t)); } @@ -3232,6 +3239,40 @@ mod geev_scaling_tests { }}; } + /// The ends of the range for a complex matrix, where faer's complex modulus overflows + /// (`2^1023 + 0i`, `2^127 + 0i`) or underflows (`2^-1074`, `2^-149`): the eigenvalues come + /// back exactly. + #[test] + fn evd_scaling_at_the_complex_range_ends() { + let top = c64::new(2.0f64.powi(1023), 0.0); + let evd = Mat::from_fn(1, 1, |_, _| top).eigen(Par::Seq).unwrap(); + assert!(evd.S()[0] == top); + let top = c32::new(2.0f32.powi(127), 0.0); + let evd = Mat::from_fn(1, 1, |_, _| top).eigen(Par::Seq).unwrap(); + assert!(evd.S()[0] == top); + // a rotation by the smallest subnormal: eigenvalues `+-i m` + let m = f64::from_bits(1); + let a = Mat::from_fn(2, 2, |i, j| match (i, j) { + (0, 1) => c64::new(-m, 0.0), + (1, 0) => c64::new(m, 0.0), + _ => c64::new(0.0, 0.0), + }); + let evd = a.eigen(Par::Seq).unwrap(); + let mut w = [evd.S()[0], evd.S()[1]]; + w.sort_by(|a, b| a.im.partial_cmp(&b.im).unwrap()); + assert!(w == [c64::new(0.0, -m), c64::new(0.0, m)]); + let m = f32::from_bits(1); + let a = Mat::from_fn(2, 2, |i, j| match (i, j) { + (0, 1) => c32::new(-m, 0.0), + (1, 0) => c32::new(m, 0.0), + _ => c32::new(0.0, 0.0), + }); + let evd = a.eigen(Par::Seq).unwrap(); + let mut w = [evd.S()[0], evd.S()[1]]; + w.sort_by(|a, b| a.im.partial_cmp(&b.im).unwrap()); + assert!(w == [c32::new(0.0, -m), c32::new(0.0, m)]); + } + #[test] fn evd_scales_its_input_into_range() { check!(f32, |z: Z| z.re as f32, |x: f32| Z::new(x as f64, 0.0), f32::EPSILON as f64, [-120, -100, 0, 100, 120], false); diff --git a/upstream-reports/faer-029-evd-input-not-scaled/README.md b/upstream-reports/faer-029-evd-input-not-scaled/README.md index 37ae9a8e..b9f93f1d 100644 --- a/upstream-reports/faer-029-evd-input-not-scaled/README.md +++ b/upstream-reports/faer-029-evd-input-not-scaled/README.md @@ -1,7 +1,7 @@ # faer-029: the nonsymmetric EVD loses every eigenvector digit when the entries are tiny -**What happens:** `Mat::eigen` returns eigenvectors whose residual `|A U - U S|` is about `1 / eps` times larger than LAPACK's, so they have no correct digit. -**Which input:** any matrix whose largest entry is below about `2^-90` (`f32`, `c32`) or `2^-960` (`f64`, `c64`). The program uses the 64 x 64 Grcar matrix scaled by a power of two. +**What happens:** `Mat::eigen` returns eigenvectors whose residual `|A U - U S|` is about `1 / eps` times larger than LAPACK's, so they have no correct digit, or, further down, eigenvectors with NaN or infinite entries, still with `Ok`. +**Which input:** a nontrivial matrix (random, Grcar and near-Jordan were tested; a diagonal one is not affected) whose largest entry is below about `2^-100` (`f32`, `c32`) or `2^-960` (`f64`, `c64`). The program uses the 64 x 64 Grcar matrix scaled by a power of two. **Which version:** `faer` 0.24.4. **Added to this list:** 2026-10-29. @@ -11,20 +11,20 @@ cargo run ``` -The program prints LAPACK's GET22 residual ratio, `|A U - U S|_1 / (|A|_1 |U|_1 eps)`, measured in `f64` on the matrix divided by its scale (exact for a power of two). +The program prints LAPACK's GET22 residual ratio (NaN when `S` or `U` has a nonfinite entry), `|A U - U S|_1 / (|A|_1 |U|_1 eps)`, measured in `f64` on the matrix divided by its scale (exact for a power of two). LAPACK's CHKHS test accepts a ratio below 20. | type | scale | expected (LAPACK `?geev`, scipy/OpenBLAS) | `faer` 0.24.4 | |---|---|---|---| | f64 | 1 | 16.4 | 10.6 | | f64 | `2^-960` | 16.4 | **1.1e15** | -| f64 | `2^-1000` | 16.4 | **9.0e14** | +| f64 | `2^-1000` | 16.4 | **`U` not finite, returned with `Ok`** | | f32 | 1 | 15.8 | 10.4 | -| f32 | `2^-100` | 15.8 | **1.8e6** | -| f32 | `2^-110` | 15.8 | **1.7e6** | +| f32 | `2^-100` | 15.8 | **`U` not finite, returned with `Ok`** | +| f32 | `2^-110` | 15.8 | **`U` not finite, returned with `Ok`** | -A sweep over `f32`, `f64`, `c32` and `c64`, n = 16 to 257, random, Grcar, near-Jordan and graded matrices, gives the same picture: from about `2^-90` (32-bit) and `2^-960` (64-bit) down, every family loses all digits, while LAPACK stays at 2 to 70. -Near the top of the range the loss is smaller: on LAPACK's own `c32` test matrices scaled to `sqrt(MAX) / 16`, `faer` gives a ratio of 42 to 59 where `cgeev` gives 29 to 40. +A sweep with this fork before its fix (which turns nonfinite eigenvectors into `NoConvergence`, so every cell has a finite ratio) over `f32`, `f64`, `c32` and `c64`, n = 16 to 257, random, Grcar and near-Jordan matrices, gives the same picture: the ratio grows from about `2^-90` (32-bit), and from `2^-100` (32-bit) and `2^-960` (64-bit) down it is `1e4` to `2e15`, while LAPACK stays at 2 to 70. +Near the top of the range the loss is smaller: on LAPACK's own `c32` test matrices scaled to `sqrt(MAX) / 16`, the fork before its fix gives a ratio of 42 to 59 where `cgeev` gives 29 to 40. ## Cause @@ -35,7 +35,7 @@ Measured stage by stage on a `c32` test matrix near the top, the Hessenberg back ## How the fork fixed it -`evd_imp` scales its input as `?GEEV` does, with a power of two instead of `smlnum / anrm`, so the scaling is exact. +`evd_imp` scales its input as `?GEEV` does, with a power of two instead of `smlnum / anrm`, so the scaling itself does not round (an entry far below the largest can still underflow when the matrix is scaled down, as with LAPACK). The largest entry is taken as the largest `|re|` or `|im|`: faer's complex modulus overflows for `2^1023 + 0i` and underflows for `2^-1074 + 0i`. The eigenvectors do not change; the eigenvalues are multiplied back. With that change this program prints 17.7 (`f64`) and 15.2 (`f32`) at every scale. diff --git a/upstream-reports/faer-029-evd-input-not-scaled/src/main.rs b/upstream-reports/faer-029-evd-input-not-scaled/src/main.rs index a88caaeb..359420e4 100644 --- a/upstream-reports/faer-029-evd-input-not-scaled/src/main.rs +++ b/upstream-reports/faer-029-evd-input-not-scaled/src/main.rs @@ -39,7 +39,9 @@ fn main() { let u = Mat::::from_fn(n, n, |i, j| c64::new(eig.U()[(i, j)].re as f64, eig.U()[(i, j)].im as f64)); (u, s, f32::EPSILON as f64) }; - let ratio = get22(&dimless, &u, &s, eps); + // `f64::max` drops a NaN, so a nonfinite result counts as a failure on its own + let finite = s.iter().all(|s| s.is_finite()) && (0..n).all(|j| (0..n).all(|i| u[(i, j)].is_finite())); + let ratio = if finite { get22(&dimless, &u, &s, eps) } else { f64::NAN }; println!("{ty:>5} | 2^{e:<4} | {:>10} | {ratio:>12.4e}", format!("< {threshold}")); failed |= !(ratio < threshold); } From 424fa12288b6cf24bf61fedf858958bdecb81889 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Thu, 29 Oct 2026 20:18:34 +0900 Subject: [PATCH 4/6] evd: the real Schur deflation-window sort always ends (#144) aggressive_early_deflation sorted the window with a bubble sort ported from LAPACK dlaqr3, but over the deflated part (dlaqr3's initial I = NS + 1 only sets the first pass's end) and with each next window set from the last swap's position. schur_swap reports success even when schur_swap_22 rejects an ill-conditioned swap, so the same pair was compared forever: Mat::eigen of the f32 matrix J_n + 1e-6 G scaled to 2^-100 did not return at n = 128 and 257. The sort now follows dlaqr3: it sorts the undeflated part infqr..ns, ends each pass at kend, one before where the previous pass stopped, and swaps with schur_swap_checked (the dlaexc ports for 1 x 2 and 2 x 1), leaving a rejected pair in place. Upstream faer 0.24.4 has the same code; the upstream report is faer-030. Closes #144 Co-Authored-By: Claude Opus 5.5 --- PATCHSET.md | 1 + faer/src/linalg/evd/schur/real_schur.rs | 181 ++++- upstream-reports/README.md | 1 + .../Cargo.lock | 761 ++++++++++++++++++ .../Cargo.toml | 13 + .../README.md | 52 ++ .../src/main.rs | 45 ++ 7 files changed, 1013 insertions(+), 41 deletions(-) create mode 100644 upstream-reports/faer-030-real-schur-aed-sort-hang/Cargo.lock create mode 100644 upstream-reports/faer-030-real-schur-aed-sort-hang/Cargo.toml create mode 100644 upstream-reports/faer-030-real-schur-aed-sort-hang/README.md create mode 100644 upstream-reports/faer-030-real-schur-aed-sort-hang/src/main.rs diff --git a/PATCHSET.md b/PATCHSET.md index 6bf94f7d..6a27da61 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -135,6 +135,7 @@ gate has changed and two kinds of adoption are now in the ledger. | Every region at the widest width that pays | The full-width gates (`par_pays`, `ParThreshold::region_pays`, `ParShape::pays_full` at the call sites, the block apply's `n_tasks`) chose between the token's whole degree and one thread, so on a wide team a region that pays at 8 or 16 ran serially. Every region now asks `utils::cost::par_width` / `ParThreshold::region_width` for the widest width that pays (every width searched, falling rows included; the GEMM kernel's `kernel_width` is the same rule) and runs on a sub-team of that width (`utils::thread::with_width`, a split into one child: no worker dispatch). Splits whose task count is capped (triangular solve jobs of 64 columns, block apply columns, `upgrade_householder_factor`'s rows per participant) search inside their cap. The factorizations that latched `par = Par::Seq` for the rest of the call (Hessenberg, tridiagonal, unblocked bidiagonal, full pivoting LU, Bunch-Kaufman, column pivoting QR) decide per step; `Never` still makes the whole call sequential; the post-loop compact-factor upgrades of the Hessenberg and tridiagonal reductions now ask their own gate (they ran under the last step's latch, i.e. sequentially). The blocked bidiagonal panel gate is a floor (sequential or the whole `par`), so the panel's products choose their own widths. The Bunch-Kaufman blocked swap loop, split over the whole degree without a gate, now has one (row swaps at the LU swap rate). `Fixed` keeps its rules (the whole degree above the floor). `par_pays`, `region_pays` and the Seq-only GEMM predicate are removed; `par-trace` records the chosen width. | EXPERIMENT (pending the maintainer's acceptance run on broccoli / evo) | tensor4all/faer#91, tensor4all/faer#95 | No timing here (performance evidence is the monitor's: broccoli 1/8/16/32/64T, evo). Tests: `utils::team::tests::{par_width_takes_the_widest_width_that_pays, with_width_runs_on_a_sub_team_of_that_width, the_families_narrow_their_regions_and_match_sequential}`, the rewritten `region_rules_follow_the_serial_time`, `fixed_and_never_keep_their_meaning`, `par_width_reads_the_degree_of_the_token`. Record: [worklog](docs/worklogs/2026-10-29-width-per-region.md). | | No heap allocation inside a parallel region (#111) | Every region on the team route allocated on the heap. `for_each_child`, behind every `for_each` and `for_each_with_par`, collected all the items into a `Vec`, one entry per column for `par_col_iter_mut`, and then a slot `Vec`. The full-pivoting LU step and both Bunch-Kaufman full-pivoting updates built an `nthreads`-long vector of pivot candidates per step, and `join_sized` boxed its two closures. The allocations were suspected of causing #109's floor-level loss (`full_piv_lu f64 128` at 8T on evo, +6.6 to +10.4%). The evo re-measurement of this change shows they were not: the cell is at 1.0025 of `main` with them gone. The change stands on its own as a cleanup of the mechanism. Now faer's dispatch and these workspaces allocate nothing. An iterator adapter that allocates in its own `split_at`, such as rayon's `map_with`, still does. `for_each_child` splits the iterator's producer on the leader into the blocks (`split_at`, remainder-relative) and each block's leader iterates its piece. Up to 64 blocks, with a producer of at most 256 bytes, are dealt from an inline array in one fork; wider or heavier splits halve recursively, one section per level. The pivot candidates come from the call's scratch: `lu_in_place_scratch` and `cholesky_in_place_scratch` grow by `par.degree()` entries. `join_sized` takes two typed slots. | KEEP (evo: no regression; `full_piv_lu f64 128` 8T at 1.0025 of `main`) | tensor4all/faer#111, tensor4all/faer#95 | Test `faer/tests/no_step_alloc.rs`: its own binary with a process-wide counting allocator. On teams of 4, 8 and 70 (the last takes the recursive path), with the regions forced parallel and a warm-up first, it counts zero allocations for full-pivoting LU 96, Bunch-Kaufman with full pivoting 96 and `join_sized`, and checks that each one dispatched. On `origin/main` the same test counts 190 allocations in one full-pivoting LU 96 call at width 4. Record: [worklog](docs/worklogs/2026-10-29-no-region-allocation.md). | | Matrix product into uninitialized storage (#136) | `matmul_with_conj_uninit(dst: MatUninitMut, lhs, conj_lhs, rhs, conj_rhs, alpha, par) -> MatMut` computes `alpha * op(lhs) * op(rhs)` into a destination of `MaybeUninit` elements and returns it as an initialized `MatMut`, for callers (strided-rs#198) that would otherwise fill the output before an overwriting product. Requested upstream as sarah-quinones/faer-rs#269; the view type `MatUninitMut` (`UninitMut`) and the first tests come from shinaoka's upstream PR #270 (same author). The destination's prior contents are never read and no `&T`/`&mut T` to an unwritten element is formed. The native types on x86_64 with AVX2 or AVX-512 go to `t4a-private-gemm-x86` with `Accum::Replace` (on a team, through `gemm_team`); other scalars write each dot product through `ptr::write`; an empty depth writes zeros. The `gemm` and `nano-gemm` crates never see uninitialized storage (`gemm`'s GEMV forms a `&mut [T]` of the destination, its Apple AMX kernel copies partial tiles of it before checking `alpha`, and `nano-gemm` scales it by `beta`), so on other targets the native types are zero-filled first and multiplied as by `matmul_with_conj`. Any strides, negative included; a batch is one view per slot. A panic leaves the storage `MaybeUninit` (possibly partly written) with nothing exposed as initialized. Existing `MatMut` APIs are unchanged. | EXPERIMENT | tensor4all/faer#136 | Tests in `linalg::matmul::tests`: NaN-prefilled destinations with padded leading dimensions for f32/f64/c32/c64, row- and column-major, conjugations, at `1x1x1` up to `130x70x90` and across the kernel's 512 depth block (`9x600x7`); the padding must stay NaN and the product must match `matmul_with_conj`. Also a team, a batched buffer, negative strides, an empty depth, a matrix-vector shape and `fx128`. Miri runs the seven small tests; six of them use `MaybeUninit::uninit()` storage (the batched one is NaN-filled), and under Miri the native types take the zero-fill route. No uninitialized read was reported. The NaN and team tests are too large for Miri. Indicative only (primerose, shared, one shape, best of 30, three runs): `2048x64x2048` f64 `Par::Seq` into existing storage is 0.89 to 0.93 of a zero fill followed by `Accum::Replace`. The saving is the output-sized fill, so it shrinks as `k` grows. Not timed on a quiet host. Record: [worklog](docs/worklogs/2026-10-29-matmul-uninit.md). | +| Real Schur: the deflation-window sort always ends (#144) | `aggressive_early_deflation` sorts the window's eigenvalues with a bubble sort ported from LAPACK `dlaqr3`, but it sorted the deflated part (`ns..jw`, from `dlaqr3`'s initial `I = NS + 1`, which only sets the first pass's end) and set each next window from the last swap's position. `schur_swap` reports success even when `schur_swap_22` rejects an ill-conditioned 2 x 2 swap, so the pass kept recording the same position and the sort never ended: `Mat::eigen` of the `f32` matrix `J_n + 1e-6 G` scaled to `2^-100` hangs at n = 128 and 257. The sort now follows `dlaqr3`: it sorts the undeflated part `infqr..ns`, ends each pass at `kend`, one before where the previous pass stopped, and swaps with `schur_swap_checked` (the `dlaexc` ports for 1 x 2 and 2 x 1), leaving a rejected pair in place. Upstream has the same code: faer 0.24.4 `src/linalg/evd/schur/real_schur.rs:921-969`. Not reported upstream. | KEEP | tensor4all/faer#144 | Correctness: the sort ends after at most `ns` passes. Test `tests::aed_sort_terminates_when_a_swap_fails` (the reproducer at n = 128 and 257 through `hessenberg_schur_real` under a 120 s timeout; it times out on `origin/main`; residual within LAPACK's `smlnum` at `2^-100` and within `10 n eps |H|` at `2^-80`, orthogonality `10 n eps`). Upstream report: [`upstream-reports/faer-030-real-schur-aed-sort-hang`](upstream-reports/faer-030-real-schur-aed-sort-hang/). | ## Ground rules diff --git a/faer/src/linalg/evd/schur/real_schur.rs b/faer/src/linalg/evd/schur/real_schur.rs index 55adc493..a2bf7936 100644 --- a/faer/src/linalg/evd/schur/real_schur.rs +++ b/faer/src/linalg/evd/schur/real_schur.rs @@ -1210,56 +1210,52 @@ fn aggressive_early_deflation( ns -= infqr; return (ns, nd); } + // Sort the undeflated eigenvalues `tw[infqr..ns]` by decreasing magnitude, + // as LAPACK `dlaqr3` does. Each pass ends at `kend`, one before where the + // previous pass stopped, so the window shrinks by at least one per pass and + // the sort ends even when a swap leaves a pair out of order. + let ev = |tw: &MatMut<'_, T>, i: usize, n: usize| -> T { + if n == 1 { + tw[(i, i)].abs() + } else { + tw[(i, i)].abs() + + tw[(i + 1, i)].abs().sqrt() * tw[(i, i + 1)].abs().sqrt() + } + }; + // Size of the block starting at `i` when the window ends at `kend`. + let block = |tw: &MatMut<'_, T>, i: usize, kend: usize| -> usize { + if i >= kend || tw[(i + 1, i)] == zero() { 1 } else { 2 } + }; let mut sorted = false; - let mut sorting_window_size = jw as isize; - while !sorted { + let mut i = ns; + while !sorted && i > infqr { sorted = true; - let mut ilst = 0isize; - let mut i1 = ns; - while i1 as isize + 1 < sorting_window_size { - let mut n1 = 1; - if tw[(i1 + 1, i1)] != zero() { - n1 = 2; - } - if i1 + n1 == jw { - ilst -= n1 as isize; - break; - } - let i2 = i1 + n1; - let mut n2 = 1; - if i2 + 1 < jw && (tw[(i2 + 1, i2)] != zero()) { - n2 = 2; - } - let (ev1, ev2); - if n1 == 1 { - ev1 = tw[(i1, i1)].abs(); - } else { - ev1 = tw[(i1, i1)].abs() - + tw[(i1 + 1, i1)].abs().sqrt() - * tw[(i1, i1 + 1)].abs().sqrt(); - } - if n2 == 1 { - ev2 = tw[(i2, i2)].abs(); - } else { - ev2 = tw[(i2, i2)].abs() - + tw[(i2 + 1, i2)].abs().sqrt() - * tw[(i2, i2 + 1)].abs().sqrt(); - } - if ev1 >= ev2 { - i1 = i2; + let kend = i - 1; + i = infqr; + let mut k = i + block(&tw, i, kend); + while k <= kend { + let n1 = k - i; + let n2 = block(&tw, k, kend); + if ev(&tw, i, n1) >= ev(&tw, k, n2) { + i = k; } else { sorted = false; - let ierr = - schur_swap(tw.rb_mut(), Some(v.rb_mut()), i1, n1, n2); + let ierr = schur_swap_checked( + tw.rb_mut(), + Some(v.rb_mut()), + i, + n1, + n2, + true, + ); if ierr == 0 { - i1 += n2; + i += n2; } else { - i1 = i2; + i = k; } - ilst = i1 as isize; } + k = i + block(&tw, i, kend); } - sorting_window_size = ilst; } let mut i = 0; while i < jw { @@ -3273,4 +3269,107 @@ mod tests { } } } + /// tensor4all/faer#144: the deflation-window sort must finish when a swap + /// is rejected. `J_n + 1e-6 G` (Jordan block with eigenvalue 1), divided by + /// its largest entry and scaled by `2^-100`, made the sort compare the same + /// pair forever in `f32` at n = 128 and 257. Each run is bounded by a + /// timeout and checked for a similarity: at `2^-100` the input is below + /// LAPACK's negligibility level `smlnum = safmin * n / eps`, so the residual + /// is held to that absolute level; the same matrix at `2^-80` is held to + /// the usual relative bound. + #[test] + fn aed_sort_terminates_when_a_swap_fails() { + use crate::linalg::evd::hessenberg::{ + hessenberg_in_place, hessenberg_in_place_scratch, + }; + use crate::linalg::evd::{ + hessenberg_schur_real, hessenberg_schur_scratch, + }; + use rand::prelude::*; + use rand_distr::StandardNormal; + use std::sync::mpsc; + use std::time::Duration; + + /// Returns `max |h - q t q^T|`, `max |h|` and `max |q^T q - I|`. + fn run(n: usize, exponent: i32) -> (f32, f32, f32) { + let rng = &mut StdRng::seed_from_u64(112); + let mut a = Mat::::from_fn(n, n, |i, j| { + let g: f32 = rng.sample(StandardNormal); + let jordan = if i == j || i + 1 == j { 1.0 } else { 0.0 }; + jordan + 1e-6 * g + }); + let scale = 2.0f32.powi(exponent) / a.norm_max(); + a *= Scale(scale); + let mut householder = Mat::::zeros(1, n - 1); + hessenberg_in_place( + a.as_mut(), + householder.as_mut(), + Par::Seq, + MemStack::new(&mut MemBuffer::new( + hessenberg_in_place_scratch::( + n, + 1, + ParShape::Seq, + Default::default(), + ), + )), + Default::default(), + ); + let h = Mat::::from_fn(n, n, |i, j| { + if i > j + 1 { 0.0 } else { a[(i, j)] } + }); + let mut t = h.clone(); + let mut q = Mat::::identity(n, n); + let mut wr = Col::::zeros(n); + let mut wi = Col::::zeros(n); + hessenberg_schur_real( + t.as_mut(), + Some(q.as_mut()), + wr.as_mut(), + wi.as_mut(), + Par::Seq, + MemStack::new(&mut MemBuffer::new( + hessenberg_schur_scratch::( + n, + true, + ParShape::Seq, + Default::default(), + ), + )), + Default::default(), + ) + .unwrap(); + assert!(t.is_all_finite() && q.is_all_finite()); + let residual = (&h - &q * &t * q.transpose()).norm_max(); + let orthogonality = + (q.transpose() * &q - Mat::::identity(n, n)).norm_max(); + (residual, h.norm_max(), orthogonality) + } + + for n in [128, 257] { + let n_eps = n as f32 * f32::EPSILON; + let smlnum = f32::MIN_POSITIVE * (n as f32 / f32::EPSILON); + for exponent in [-100, -80] { + let (send, recv) = mpsc::channel(); + std::thread::spawn(move || { + send.send(run(n, exponent)).unwrap() + }); + let (residual, norm, orthogonality) = recv + .recv_timeout(Duration::from_secs(120)) + .unwrap_or_else(|_| { + panic!("n = {n}, 2^{exponent}: real Schur did not finish") + }); + let bound = if exponent == -100 { + smlnum + } else { + 10.0 * n_eps * norm + }; + assert!(residual <= bound, "n = {n}, 2^{exponent}: residual {residual}, bound {bound}"); + assert!( + orthogonality <= 10.0 * n_eps, + "n = {n}, 2^{exponent}: orthogonality {orthogonality}" + ); + } + } + } } diff --git a/upstream-reports/README.md b/upstream-reports/README.md index 7224376c..e4fcd4ac 100644 --- a/upstream-reports/README.md +++ b/upstream-reports/README.md @@ -43,6 +43,7 @@ The program prints the expected and the actual value, and exits with status 1 wh | [faer-027](faer-027-qr-near-dependent-column-skip/) | Unpivoted QR skips a nearly dependent column of a tall matrix | `faer` 0.24.4 | inaccurate value | tensor4all/faer#1 (`2e6ea52`, `d16929b`) | not yet | 2026-10-29 | | [faer-028](faer-028-sparse-lu-numeric-singular/) | Sparse LU factors a numerically singular matrix without an error | `faer` 0.24.4 | NaN | tensor4all/faer#132 | not yet | 2026-10-29 | | [faer-029](faer-029-evd-input-not-scaled/) | The nonsymmetric EVD loses every eigenvector digit when the entries are tiny | `faer` 0.24.4 | inaccurate value | tensor4all/faer#112 | not yet | 2026-10-29 | +| [faer-030](faer-030-real-schur-aed-sort-hang/) | Real Schur QR does not return when the deflation-window sort cannot swap a pair | `faer` 0.24.4 | hang | tensor4all/faer#144 | not yet | 2026-10-29 | | [gemm-x86-001](gemm-x86-001-diag-negative-stride/) | A product with a diagonal gives a wrong value when the left matrix is stored backwards | `private-gemm-x86` 0.1.22 | wrong value | tensor4all/t4a-private-gemm-x86#5 (`84ca6db`) | not yet | 2026-10-29 | "Fixed upstream" was checked against upstream `main`: `faer` at `0ecd846` (2026-10-22) and `private-gemm-x86` at `71872b9` (its 0.1.22 release). diff --git a/upstream-reports/faer-030-real-schur-aed-sort-hang/Cargo.lock b/upstream-reports/faer-030-real-schur-aed-sort-hang/Cargo.lock new file mode 100644 index 00000000..b3aa3b90 --- /dev/null +++ b/upstream-reports/faer-030-real-schur-aed-sort-hang/Cargo.lock @@ -0,0 +1,761 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 4 + +[[package]] +name = "autocfg" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" + +[[package]] +name = "bitflags" +version = "2.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3ded4057c258ba199e2d26386d3af3780957ecaee6c4ef4041c6b4b8b97c0b06" + +[[package]] +name = "bytemuck" +version = "1.25.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" +dependencies = [ + "bytemuck_derive", +] + +[[package]] +name = "bytemuck_derive" +version = "1.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6a1f896587b6f2c069c73d2f0913e2d590c3990285cd2f0b6aa02b786b4c679c" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.9", +] + +[[package]] +name = "byteorder" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" + +[[package]] +name = "cfg-if" +version = "1.0.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4e7648175b45a9a48536d676f68d918270699102aa8dab5496df06904c914600" + +[[package]] +name = "crunchy" +version = "0.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "460fbee9c2c2f33933d720630a6a0bac33ba7053db5344fac858d4b8952d77d5" + +[[package]] +name = "defer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "930c7171c8df9fb1782bdf9b918ed9ed2d33d1d22300abb754f9085bc48bf8e8" + +[[package]] +name = "dyn-stack" +version = "0.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1c4713e43e2886ba72b8271aa66c93d722116acf7a75555cce11dcde84388fe8" +dependencies = [ + "bytemuck", + "dyn-stack-macros", +] + +[[package]] +name = "dyn-stack-macros" +version = "0.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e1d926b4d407d372f141f93bb444696142c29d32962ccbd3531117cf3aa0bfa9" + +[[package]] +name = "enum-as-inner" +version = "0.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a1e6a265c649f3f5979b601d26f1d05ada116434c87741c9493cb56218f76cbc" +dependencies = [ + "heck", + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "equator" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c35da53b5a021d2484a7cc49b2ac7f2d840f8236a286f84202369bd338d761ea" +dependencies = [ + "equator-macro 0.2.1", +] + +[[package]] +name = "equator" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "02da895aab06bbebefb6b2595f6d637b18c9ff629b4cd840965bb3164e4194b0" +dependencies = [ + "equator-macro 0.6.0", +] + +[[package]] +name = "equator-macro" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3bf679796c0322556351f287a51b49e48f7c4986e727b5dd78c972d30e2e16cc" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "equator-macro" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b14b339eb76d07f052cdbad76ca7c1310e56173a138095d3bf42a23c06ef5d8" + +[[package]] +name = "faer" +version = "0.24.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5ab6df3dd147fe8d702a288b95bcd8fcc499ab572fc80da6828f60cd4d524d67" +dependencies = [ + "bytemuck", + "dyn-stack", + "equator 0.6.0", + "faer-traits", + "gemm", + "generativity", + "libm", + "nano-gemm", + "num-complex", + "num-traits", + "private-gemm-x86", + "pulp", + "reborrow", +] + +[[package]] +name = "faer-030-real-schur-aed-sort-hang" +version = "0.0.0" +dependencies = [ + "faer", + "rand", + "rand_distr", +] + +[[package]] +name = "faer-traits" +version = "0.24.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b87d23ed7ab1f26c0cba0e5b9e061a796fbb7dc170fa8bee6970055a1308bb0f" +dependencies = [ + "bytemuck", + "dyn-stack", + "generativity", + "libm", + "num-complex", + "num-traits", + "pulp", + "qd", + "reborrow", +] + +[[package]] +name = "gemm" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa0673db364b12263d103b68337a68fbecc541d6f6b61ba72fe438654709eacb" +dependencies = [ + "dyn-stack", + "gemm-c32", + "gemm-c64", + "gemm-common", + "gemm-f16", + "gemm-f32", + "gemm-f64", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-c32" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "086936dbdcb99e37aad81d320f98f670e53c1e55a98bee70573e83f95beb128c" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-c64" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "20c8aeeeec425959bda4d9827664029ba1501a90a0d1e6228e48bef741db3a3f" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-common" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "88027625910cc9b1085aaaa1c4bc46bb3a36aad323452b33c25b5e4e7c8e2a3e" +dependencies = [ + "bytemuck", + "dyn-stack", + "half", + "libm", + "num-complex", + "num-traits", + "once_cell", + "paste", + "pulp", + "raw-cpuid", + "seq-macro", + "sysctl", +] + +[[package]] +name = "gemm-f16" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e3df7a55202e6cd6739d82ae3399c8e0c7e1402859b30e4cb780e61525d9486e" +dependencies = [ + "dyn-stack", + "gemm-common", + "gemm-f32", + "half", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-f32" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "02e0b8c9da1fbec6e3e3ab2ce6bc259ef18eb5f6f0d3e4edf54b75f9fd41a81c" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-f64" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "056131e8f2a521bfab322f804ccd652520c79700d81209e9d9275bbdecaadc6a" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "generativity" +version = "1.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d2c81fb5260e37854d09d5c87183309fd8c555b75289427884b25660bc87a85e" + +[[package]] +name = "getrandom" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "899def5c37c4fd7b2664648c28120ecec138e4d395b459e5ca34f9cce2dd77fd" +dependencies = [ + "cfg-if", + "libc", + "r-efi", + "wasip2", +] + +[[package]] +name = "half" +version = "2.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6ea2d84b969582b4b1864a92dc5d27cd2b77b622a8d79306834f1be5ba20d84b" +dependencies = [ + "bytemuck", + "cfg-if", + "crunchy", + "num-traits", + "zerocopy", +] + +[[package]] +name = "heck" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" + +[[package]] +name = "hermit-abi" +version = "0.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e17592d60ebacc7d5e169f4663c5f84f9161cc90328abcfe8456f41e4dfcb284" + +[[package]] +name = "interpol" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eb58032ba748f4010d15912a1855a8a0b1ba9eaad3395b0c171c09b3b356ae50" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "libc" +version = "0.2.191" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "65cd87c9a6d60536580b316ad985726104730a5151c99845f2490428ef91240a" + +[[package]] +name = "libm" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6d2cec3eae94f9f509c767b45932f1ada8350c4bdb85af2fcab4a3c14807981" + +[[package]] +name = "nano-gemm" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9e04345dc84b498ff89fe0d38543d1f170da9e43a2c2bcee73a0f9069f72d081" +dependencies = [ + "equator 0.2.2", + "nano-gemm-c32", + "nano-gemm-c64", + "nano-gemm-codegen", + "nano-gemm-core", + "nano-gemm-f32", + "nano-gemm-f64", + "num-complex", +] + +[[package]] +name = "nano-gemm-c32" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0775b1e2520e64deee8fc78b7732e3091fb7585017c0b0f9f4b451757bbbc562" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", + "num-complex", +] + +[[package]] +name = "nano-gemm-c64" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9af49a20d58816e6b5ee65f64142e50edb5eba152678d4bb7377fcbf63f8437a" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", + "num-complex", +] + +[[package]] +name = "nano-gemm-codegen" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6cc8d495c791627779477a2cf5df60049f5b165342610eb0d76bee5ff5c5d74c" + +[[package]] +name = "nano-gemm-core" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d998dfa644de87a0f8660e5ea511d7cb5c33b5a2d9847b7af57a2565105089f0" + +[[package]] +name = "nano-gemm-f32" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "879d962e79bc8952e4ad21ca4845a21132540ed3f5e01184b2ff7f720e666523" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", +] + +[[package]] +name = "nano-gemm-f64" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9a513473dce7dc00c7e7c318481ca4494034e76997218d8dad51bd9f007a815" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", +] + +[[package]] +name = "num-complex" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" +dependencies = [ + "bytemuck", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", + "libm", +] + +[[package]] +name = "num_cpus" +version = "1.17.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91df4bbde75afed763b708b7eee1e8e7651e02d97f6d5dd763e89367e957b23b" +dependencies = [ + "hermit-abi", + "libc", +] + +[[package]] +name = "once_cell" +version = "1.21.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50" + +[[package]] +name = "paste" +version = "1.0.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" + +[[package]] +name = "ppv-lite86" +version = "0.2.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85eae3c4ed2f50dcfe72643da4befc30deadb458a9b590d720cde2f2b1e97da9" +dependencies = [ + "zerocopy", +] + +[[package]] +name = "private-gemm-x86" +version = "0.1.22" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bba0a7831d1297712e7e7b33759b0c837941e79196dd05df506d5f36f7f78619" +dependencies = [ + "defer", + "interpol", + "num_cpus", + "raw-cpuid", +] + +[[package]] +name = "proc-macro2" +version = "1.0.107" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "985e7ec9bb745e6ce6535b544d84d6cd6f7ad8bd711c398938ae983b91a766d9" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "pulp" +version = "0.22.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "046aa45b989642ec2e4717c8e72d677b13edd831a4d3b6cf37d9a3e54912496a" +dependencies = [ + "bytemuck", + "cfg-if", + "libm", + "num-complex", + "paste", + "pulp-wasm-simd-flag", + "raw-cpuid", + "reborrow", + "version_check", +] + +[[package]] +name = "pulp-wasm-simd-flag" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1d8f70e07b9c3962945a74e59ca1c511bba65b6419468acc217c457d93f3c740" + +[[package]] +name = "qd" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "15f1304a5aecdcfe9ee72fbba90aa37b3aa067a69d14cb7f3d9deada0be7c07c" +dependencies = [ + "bytemuck", + "libm", + "num-traits", + "pulp", +] + +[[package]] +name = "quote" +version = "1.0.47" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fbf4db142a473a8d80c26bbf18454ed458bf8d26c8219c331daecfdbd079001" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "r-efi" +version = "5.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "69cdb34c158ceb288df11e18b4bd39de994f6657d83847bdffdbd7f346754b0f" + +[[package]] +name = "rand" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9ef1d0d795eb7d84685bca4f72f3649f064e6641543d3a8c415898726a57b41" +dependencies = [ + "rand_chacha", + "rand_core", +] + +[[package]] +name = "rand_chacha" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d3022b5f1df60f26e1ffddd6c66e8aa15de382ae63b3a0c1bfc0e4d3e3f325cb" +dependencies = [ + "ppv-lite86", + "rand_core", +] + +[[package]] +name = "rand_core" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "76afc826de14238e6e8c374ddcc1fa19e374fd8dd986b0d2af0d02377261d83c" +dependencies = [ + "getrandom", +] + +[[package]] +name = "rand_distr" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6a8615d50dcf34fa31f7ab52692afec947c4dd0ab803cc87cb3b0b4570ff7463" +dependencies = [ + "num-traits", + "rand", +] + +[[package]] +name = "raw-cpuid" +version = "11.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "498cd0dc59d73224351ee52a95fee0f1a617a2eae0e7d9d720cc622c73a54186" +dependencies = [ + "bitflags", +] + +[[package]] +name = "reborrow" +version = "0.5.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "03251193000f4bd3b042892be858ee50e8b3719f2b08e5833ac4353724632430" + +[[package]] +name = "same-file" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93fc1dc3aaa9bfed95e02e6eadabb4baf7e3078b0bd1b4d7b6b0b68378900502" +dependencies = [ + "winapi-util", +] + +[[package]] +name = "seq-macro" +version = "0.3.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1bc711410fbe7399f390ca1c3b60ad0f53f80e95c5eb935e52268a0e2cd49acc" + +[[package]] +name = "syn" +version = "1.0.109" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "2.0.119" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "872831b642d1a07999a962a351ed35b955ea2cfc8f3862091e2a240a84f17297" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "3.0.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d78c8dee4c7bf0e14673097256fed6142ce9d3b85a408189d07482442145823b" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "sysctl" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "01198a2debb237c62b6826ec7081082d951f46dbb64b0e8c7649a452230d1dfc" +dependencies = [ + "bitflags", + "byteorder", + "enum-as-inner", + "libc", + "thiserror", + "walkdir", +] + +[[package]] +name = "thiserror" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6aaf5339b578ea85b50e080feb250a3e8ae8cfcdff9a461c9ec2904bc923f52" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "unicode-ident" +version = "1.0.27" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a2c754d6c33795a1c324727428e5a7dedb5b06195f9890bdbcba760d3e246563" + +[[package]] +name = "version_check" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" + +[[package]] +name = "walkdir" +version = "2.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "29790946404f91d9c5d06f9874efddea1dc06c5efe94541a7d6863108e3a5e4b" +dependencies = [ + "same-file", + "winapi-util", +] + +[[package]] +name = "wasip2" +version = "1.0.4+wasi-0.2.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b67efb37e106e55ce722a510d6b5f9c17f083e5fc79afc2badeb12cc313d9487" +dependencies = [ + "wit-bindgen", +] + +[[package]] +name = "winapi-util" +version = "0.1.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c2a7b1c03c876122aa43f3020e6c3c3ee5c05081c9a00739faf7503aeba10d22" +dependencies = [ + "windows-sys", +] + +[[package]] +name = "windows-link" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f0805222e57f7521d6a62e36fa9163bc891acd422f971defe97d64e70d0a4fe5" + +[[package]] +name = "windows-sys" +version = "0.61.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ae137229bcbd6cdf0f7b80a31df61766145077ddf49416a728b02cb3921ff3fc" +dependencies = [ + "windows-link", +] + +[[package]] +name = "wit-bindgen" +version = "0.57.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ebf944e87a7c253233ad6766e082e3cd714b5d03812acc24c318f549614536e" + +[[package]] +name = "zerocopy" +version = "0.8.63" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e5fe1f8f1b06191a00962174c61aa5005e0bb391a6d80d07e24d115c01a92ed8" +dependencies = [ + "zerocopy-derive", +] + +[[package]] +name = "zerocopy-derive" +version = "0.8.63" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "863ad3ac83293fb4d740aedbfdc9240dd8d1a50c1099acd76ce80ce7c7230c7f" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] diff --git a/upstream-reports/faer-030-real-schur-aed-sort-hang/Cargo.toml b/upstream-reports/faer-030-real-schur-aed-sort-hang/Cargo.toml new file mode 100644 index 00000000..ab29dc54 --- /dev/null +++ b/upstream-reports/faer-030-real-schur-aed-sort-hang/Cargo.toml @@ -0,0 +1,13 @@ +[package] +name = "faer-030-real-schur-aed-sort-hang" +version = "0.0.0" +edition = "2021" +publish = false + +# Not part of the t4a-faer workspace: this reproduces the published crate on its own. +[workspace] + +[dependencies] +faer = { version = "=0.24.4", default-features = false, features = ["std", "linalg"] } +rand = "0.9" +rand_distr = "0.5" diff --git a/upstream-reports/faer-030-real-schur-aed-sort-hang/README.md b/upstream-reports/faer-030-real-schur-aed-sort-hang/README.md new file mode 100644 index 00000000..a9a0bdbe --- /dev/null +++ b/upstream-reports/faer-030-real-schur-aed-sort-hang/README.md @@ -0,0 +1,52 @@ +# faer-030: real Schur QR does not return when the deflation-window sort cannot swap a pair + +**What happens:** `Mat::eigen` never returns for some finite `f32` matrices; the thread spins in the real Schur QR. +**Which input:** the near-Jordan matrix `J_n + 1e-6 G` (Jordan block with eigenvalue 1, `G` standard normal from `StdRng` seed 112), divided by its largest entry and scaled by `2^-100`, at n = 128 and n = 257. +**Which version:** `faer` 0.24.4. +**Added to this list:** 2026-10-29. + +## Reproduce + +```sh +cargo run --release +``` + +The program runs `Mat::eigen` on that matrix for n = 128 and n = 257 in a separate thread and waits 60 s for each. + +| n | expected | `faer` 0.24.4 | +|---|---|---| +| 128 | returns | **no return after 60 s** | +| 257 | returns | **no return after 60 s** | + +The same matrices at scales `2^0` to `2^-90` return. +A backtrace taken while it spins ends in `schur_swap_22` called from `aggressive_early_deflation`. + +## Cause + +`aggressive_early_deflation` sorts the eigenvalues of the deflation window with a bubble sort (`src/linalg/evd/schur/real_schur.rs:921-969`). +It is a port of LAPACK `dlaqr3`, with two differences: + +* `dlaqr3` sorts the undeflated part `T(INFQR+1:NS)`. The port sorts from `ns` to `jw`, the deflated part: it took `dlaqr3`'s initial `I = NS + 1`, which only sets the first pass's end `KEND = I - 1`, as the start of the sort. +* `dlaqr3` ends each pass at `KEND`, one before where the previous pass stopped, so the window shrinks by at least one per pass and the sort ends whatever the swaps do. The port sets the next window from the position of the last swap. + +`schur_swap_22` rejects an ill-conditioned swap of two 2 x 2 blocks and leaves the matrix unchanged, but `schur_swap` reports success either way. +The pass then records the same swap position as the next window, the pair is still out of order, and the next pass makes the same attempt, forever. + +## How the fork fixed it + +The sort follows `dlaqr3`: it sorts the undeflated part and bounds each pass by `KEND`. +It swaps with the checked swap, which uses the `dlaexc` ports for the 1 x 2 and 2 x 1 cases, and a rejected swap leaves the pair in place, as `dlaqr3` does after `DTREXC` fails. +With that change the fork returns for both sizes within a few seconds. +At this scale the input lies below LAPACK's negligibility level `smlnum = safmin * n / eps`, so the Schur form reproduces it only to within `smlnum` in absolute terms, as LAPACK's `?HSEQR` does for unscaled input; the same matrix scaled to `2^-80` gives the usual small relative residual. + +Fork regression test, in `faer/src/linalg/evd/schur/real_schur.rs`: `tests::aed_sort_terminates_when_a_swap_fails`. + +## Affected versions + +Confirmed by running this program: `faer` 0.24.4. +Not fixed upstream: the same code is in upstream `main` at `0ecd846` (2026-10-22). + +## References + +* Fork fix: tensor4all/faer#144 (PR to follow). +* Where it was found: tensor4all/faer#112, the LAPACK comparison sweep. diff --git a/upstream-reports/faer-030-real-schur-aed-sort-hang/src/main.rs b/upstream-reports/faer-030-real-schur-aed-sort-hang/src/main.rs new file mode 100644 index 00000000..04714d68 --- /dev/null +++ b/upstream-reports/faer-030-real-schur-aed-sort-hang/src/main.rs @@ -0,0 +1,45 @@ +// `Mat::eigen` of an f32 near-Jordan matrix, J_n + 1e-6 G (Jordan block with eigenvalue 1, +// G standard normal from StdRng seed 112), divided by its largest entry and scaled by 2^-100. +// All entries are finite, normal f32 values. +// The call should return in well under a second; the program gives each size 60 s. +use faer::Mat; +use rand::prelude::*; +use rand_distr::StandardNormal; +use std::sync::mpsc; +use std::time::{Duration, Instant}; + +fn input(n: usize) -> Mat { + let rng = &mut StdRng::seed_from_u64(112); + let a = Mat::::from_fn(n, n, |i, j| { + let g: f32 = rng.sample(StandardNormal); + let jordan = if i == j || i + 1 == j { 1.0 } else { 0.0 }; + jordan + 1e-6 * g + }); + let max = a.norm_max(); + Mat::from_fn(n, n, |i, j| a[(i, j)] / max * 2f32.powi(-100)) +} + +fn main() { + let mut bad = false; + for n in [128, 257] { + let a = input(n); + let (send, recv) = mpsc::channel(); + let start = Instant::now(); + std::thread::spawn(move || { + let ok = a.eigen().is_ok(); + let _ = send.send(ok); + }); + match recv.recv_timeout(Duration::from_secs(60)) { + Ok(ok) => println!("n = {n}: returned after {:.2} s (Ok: {ok})", start.elapsed().as_secs_f64()), + Err(_) => { + println!("n = {n}: no return after 60 s"); + bad = true; + } + } + } + println!("expected: both calls return"); + if bad { + println!("WRONG: Mat::eigen does not return"); + std::process::exit(1); + } +} From f614ccbea7b132f7a6bb9a578406e816641aba1f Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Thu, 29 Oct 2026 20:08:18 +0900 Subject: [PATCH 5/6] evd: the 2 x 2 block solves of the real eigenvector back substitution follow ?LALN2 (#145) The back substitution for a complex eigenvalue of a real Schur form solved each 2 x 2 diagonal block through its inverse and compared the block's determinant with the entry-sized clamp `eps_norm`. A determinant is a product of two entries, so every block with entries below `sqrt(eps_norm)` got the determinant `eps_norm` with its original numerators. After `scaled_into_range` the blocks of a strongly graded matrix are all that small, so every complex pair got a wrong eigenvector (GET22 up to 2.6e9). Stage comparison with LAPACK 3.12.0 on f64 n = 64: the Schur form is fine (backward error 8.0); on that same Schur form `dtrevc3` gives 2.2e-4 and faer 3.4e7, and every failing eigen-block is a complex pair. Balancing is not the cause (`dgeev` without it: 0.58; fixed faer: 0.45). Both 2 x 2 solves now use `?LALN2`'s method (`solve_shifted_2x2_cplx`): `eps_norm * I` when every entry is below `eps_norm`, otherwise elimination with complete pivoting and the second pivot raised to `eps_norm`. That covers the real-shift solve too, whose inverse through `a^2 - b c` underflowed to `inf` for tiny blocks. Raising only the determinant was rejected in review: with the original numerators it is not a solved system (4.5e14 on a 4 x 4 Schur form, now a test). Upstream faer 0.24.4 has the same comparison; it is added to upstream-reports as faer-031 with a 4 x 4 input that fails there and passes with this solve. Closes #145 Co-Authored-By: Claude Opus 5.5 --- PATCHSET.md | 1 + faer/src/linalg/evd/mod.rs | 122 +++- .../evd/tests/lapack_nonsymmetric_evd.rs | 131 ++++ upstream-reports/README.md | 1 + .../faer-031-real-evd-pair-clamp/Cargo.lock | 678 ++++++++++++++++++ .../faer-031-real-evd-pair-clamp/Cargo.toml | 11 + .../faer-031-real-evd-pair-clamp/README.md | 56 ++ .../faer-031-real-evd-pair-clamp/src/main.rs | 39 + 8 files changed, 1012 insertions(+), 27 deletions(-) create mode 100644 upstream-reports/faer-031-real-evd-pair-clamp/Cargo.lock create mode 100644 upstream-reports/faer-031-real-evd-pair-clamp/Cargo.toml create mode 100644 upstream-reports/faer-031-real-evd-pair-clamp/README.md create mode 100644 upstream-reports/faer-031-real-evd-pair-clamp/src/main.rs diff --git a/PATCHSET.md b/PATCHSET.md index 6a27da61..fedee0bc 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -102,6 +102,7 @@ gate has changed and two kinds of adoption are now in the ledger. | 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/). | +| The 2 x 2 block solves of the real eigenvector back substitution follow `?LALN2` (#145) | `solve_cplx_shifted_upper_quasi_triangular_system` solved each 2 x 2 diagonal block of the real Schur form for a complex eigenvalue through the block's inverse, and compared its determinant with the entry-sized clamp `eps_norm`. The determinant has the size of a product of two entries, so every block with entries below `sqrt(eps_norm)` got the determinant `eps_norm` with its original numerators. After `scaled_into_range` (#112) the blocks of a strongly graded matrix `D A D^-1` (`D = diag(10^linspace(-6, 6, n))`) are all that small, and every complex pair got a wrong eigenvector: GET22 `1.0e9` at `f64` n = 64 in the issue's sweep. Stage comparison against LAPACK 3.12.0 (`f64` n = 64): the reduction is fine (Schur backward error 8.0 eps |A|); on faer's own Schur form `dtrevc3` gives `2.2e-4` and faer `3.4e7`, and all 24 failing eigen-blocks are complex pairs. Both 2 x 2 solves (complex shift, and the real-shift one above a real eigenvalue, which inverted through `a^2 - b c` and gave `inf` when that underflowed) now use `solve_shifted_2x2_cplx`, `?LALN2`'s method: `eps_norm * I` when every entry is below `eps_norm`, otherwise elimination with complete pivoting and the second pivot raised to `eps_norm`. Raising only the determinant was tried first and rejected in review: with the original numerators it is not a solved system (residual `4.5e14` on a 4 x 4 Schur form with a small determinant next to large entries). Balancing (`?GEBAL`) is not the cause: `dgeev` without balancing gives 0.58 on the same matrix, the fixed faer 0.45; balancing would bring both to about `1e-10`. Upstream has the same comparison (`z < eps * norm`, faer 0.24.4 `src/linalg/evd/mod.rs:828`); without the fork's scaling it fails on smaller inputs, for example two small, close complex pairs coupled by large entries. | KEEP | tensor4all/faer#145 | Correctness, no timing claim: one 2 x 2 elimination per block in place of one inverse. Tests (four dtypes, the real ones reach the changed code): `lapack_nonsymmetric_evd_graded_pairs_keep_their_residual` (graded n = 16, 64, 128 and the coupled 4 x 4, left and right, GET22 against `get22_threshold`; before: `5.9e8` at n = 16, `1.5e12` on the 4 x 4), `lapack_nonsymmetric_evd_pair_solve_with_a_small_determinant` (fails with the determinant floor: `4.5e14`; on `main`: `f32` `4.2e6`) and `lapack_nonsymmetric_evd_real_shift_block_below_the_square_root_of_min_positive` (`NoConvergence` on `main`). `linalg::evd` with `--include-ignored`: the one failure is the known c32 kinds 6/7 case of #112. Upstream report: [`upstream-reports/faer-031-real-evd-pair-clamp`](upstream-reports/faer-031-real-evd-pair-clamp/). | | 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). | diff --git a/faer/src/linalg/evd/mod.rs b/faer/src/linalg/evd/mod.rs index 8347365f..6b85d6d0 100644 --- a/faer/src/linalg/evd/mod.rs +++ b/faer/src/linalg/evd/mod.rs @@ -901,14 +901,19 @@ fn solve_real_shifted_upper_quasi_triangular_system( ); x[i - 1] -= dot0; x[i] -= dot1; - let ref a = &A[(i, i)] - p; - let ref b = A[(i - 1, i)]; - let ref c = A[(i, i - 1)]; - let ref r0 = x[i - 1]; - let ref r1 = x[i]; - let ref inv_det = (a.abs2() - b * c).recip(); - (x[i - 1], x[i]) = - ((a * r0 - b * r1) * inv_det, (a * r1 - c * r0) * inv_det); + // The same `?LALN2` solve as for a complex shift, with a zero + // imaginary part. The inverse through `a^2 - b c` underflowed to a + // zero determinant for a block whose entries are below the square + // root of the smallest positive number, and returned `inf`. + let ((y0, _), (y1, _)) = solve_shifted_2x2_cplx( + (&A[(i, i)] - p, zero()), + A[(i - 1, i)].copy(), + A[(i, i - 1)].copy(), + (x[i - 1].copy(), zero()), + (x[i].copy(), zero()), + &eps_norm, + ); + (x[i - 1], x[i]) = (y0, y1); let a0 = x[i - 1].abs(); let a1 = x[i].abs(); let m = if a0 > a1 { a0 } else { a1 }; @@ -966,6 +971,70 @@ fn solve_real_shifted_upper_quasi_triangular_system( f0 * f1 } } +/// Solves `[[a, b], [c, a]] [x0, x1]^T = [r0, r1]^T` for complex `a` (as `(re, im)`), real +/// `b`, `c` and complex right-hand sides, as LAPACK's `?LALN2` does for a complex shift of a +/// 2 x 2 block: when every entry is below `smin` in `|re| + |im|`, the block is taken as +/// `smin * I`; otherwise Gaussian elimination with complete pivoting, with the second pivot +/// raised to `smin` when it is smaller. The caller rescales a large solution. +fn solve_shifted_2x2_cplx( + a: (T, T), + b: T, + c: T, + r0: (T, T), + r1: (T, T), + smin: &T, +) -> ((T, T), (T, T)) { + type C = (T, T); + fn abs1(z: &C) -> T { + z.0.abs() + z.1.abs() + } + fn mul(x: &C, y: &C) -> C { + (&x.0 * &y.0 - &x.1 * &y.1, &x.0 * &y.1 + &x.1 * &y.0) + } + fn sub(x: &C, y: &C) -> C { + (&x.0 - &y.0, &x.1 - &y.1) + } + fn div(x: &C, y: &C) -> C { + // scaled by the larger component of `y`, so `|y|^2` neither overflows nor underflows + let s = y.0.abs().fmax(y.1.abs()); + let (yr, yi) = (&y.0 / &s, &y.1 / &s); + let (xr, xi) = (&x.0 / &s, &x.1 / &s); + let d = &yr * &yr + &yi * &yi; + ((&xr * &yr + &xi * &yi) / &d, (&xi * &yr - &xr * &yi) / &d) + } + let zero = zero::; + let m = [[(a.0.copy(), a.1.copy()), (b, zero())], [(c, zero()), a]]; + let (mut p, mut q) = (0usize, 0usize); + let mut cmax = abs1(&m[0][0]); + for (i, j) in [(0, 1), (1, 0)] { + let v = abs1(&m[i][j]); + if v > cmax { + cmax = v; + (p, q) = (i, j); + } + } + if cmax < *smin { + let inv = smin.recip(); + return ( + (&r0.0 * &inv, &r0.1 * &inv), + (&r1.0 * &inv, &r1.1 * &inv), + ); + } + // pivot row `p`, pivot column `q`; the other row and column are `1 - p`, `1 - q` + let (rp, ro) = if p == 0 { (r0, r1) } else { (r1, r0) }; + let u11 = (m[p][q].0.copy(), m[p][q].1.copy()); + let u12 = (m[p][1 - q].0.copy(), m[p][1 - q].1.copy()); + let l21 = div(&m[1 - p][q], &u11); + let mut u22 = sub(&m[1 - p][1 - q], &mul(&l21, &u12)); + if abs1(&u22) < *smin { + u22 = (smin.copy(), zero()); + } + let y2 = sub(&ro, &mul(&l21, &rp)); + let x2 = div(&y2, &u22); + let x1 = div(&sub(&rp, &mul(&u12, &x2)), &u11); + // unknown `q` is `x1`, the other is `x2` + if q == 0 { (x1, x2) } else { (x2, x1) } +} fn solve_cplx_shifted_upper_quasi_triangular_system( A: MatRef<'_, T>, p: T, @@ -1043,25 +1112,24 @@ fn solve_cplx_shifted_upper_quasi_triangular_system( let ref r0_im = x[(i - 1, 1)]; let ref r1_re = x[(i, 0)]; let ref r1_im = x[(i, 1)]; - let mut z_re = a_re.abs2() - a_im.abs2() - b * c; - let mut z_im = (a_re * a_im).mul_pow2(from_f64::(2.0)); - let mut z = z_re.hypot(&z_im); - if z < eps_norm { - z_re = eps_norm.copy(); - z_im = zero(); - z = z_re.copy(); - } - let ref z = z; - let ref z_re = (z_re / z) / z; - let ref z_im = (-z_im / z) / z; - let ref x0_re = (a_re * r0_re - a_im * r0_im) - b * r1_re; - let ref x0_im = (a_re * r0_im + a_im * r0_re) - b * r1_im; - let ref x1_re = (a_re * r1_re - a_im * r1_im) - c * r0_re; - let ref x1_im = (a_re * r1_im + a_im * r1_re) - c * r0_im; - x[(i - 1, 0)] = x0_re * z_re - x0_im * z_im; - x[(i - 1, 1)] = x0_re * z_im + x0_im * z_re; - x[(i, 0)] = x1_re * z_re - x1_im * z_im; - x[(i, 1)] = x1_re * z_im + x1_im * z_re; + // `?LALN2`'s complex 2 x 2 solve (#145): complete pivoting, the + // second pivot kept at least `eps_norm`, and `eps_norm * I` for a + // block whose entries are all below it. Inverting the block through + // its determinant compared that determinant (a product of two + // entries) with the entry-sized `eps_norm`, and clamping only the + // determinant does not give a perturbed system that is solved. + let ((y0_re, y0_im), (y1_re, y1_im)) = solve_shifted_2x2_cplx( + (a_re.copy(), a_im.copy()), + b.copy(), + c.copy(), + (r0_re.copy(), r0_im.copy()), + (r1_re.copy(), r1_im.copy()), + &eps_norm, + ); + x[(i - 1, 0)] = y0_re; + x[(i - 1, 1)] = y0_im; + x[(i, 0)] = y1_re; + x[(i, 1)] = y1_im; let mut m = x[(i - 1, 0)].abs(); let v = x[(i - 1, 1)].abs(); if v > m { diff --git a/faer/src/linalg/evd/tests/lapack_nonsymmetric_evd.rs b/faer/src/linalg/evd/tests/lapack_nonsymmetric_evd.rs index f0d44b80..fbee70b6 100644 --- a/faer/src/linalg/evd/tests/lapack_nonsymmetric_evd.rs +++ b/faer/src/linalg/evd/tests/lapack_nonsymmetric_evd.rs @@ -577,3 +577,134 @@ fn lapack_nonsymmetric_evd_clustered_residuals_do_not_grow_with_n() { } } } +/// A strongly graded matrix `D A D^-1`, `D = diag(10^linspace(-6, 6, n))` (#145). Its Schur form +/// spans many orders of magnitude, so after `scaled_into_range` the 2 x 2 blocks of the +/// complex pairs have entries far below one. The back substitution compared the determinant of +/// such a block with the entry-sized clamp, replaced the determinant of every block and returned +/// a wrong vector for each complex pair: a GET22 ratio of `5.9e8` here at `f64` n = 16. On the `f64` n = 64 Schur form, +/// LAPACK's `dtrevc3` gives `2.2e-4` and the fixed back substitution `3.1e-4`. +/// +/// The second matrix also fails upstream: two complex pairs `+-1e-5 i` and `+-2e-5 i` coupled +/// by `1e8`, so the shifted block's determinant (about `1e-10`) is below an entry-sized clamp +/// although its entries are not. Before the fix: `1.5e12` (`f64`); `dgeev` gives `9.6e-14`. +#[test] +fn lapack_nonsymmetric_evd_graded_pairs_keep_their_residual() { + let coupled = Mat::from_fn(4, 4, |i, j| { + scalar( + match (i, j) { + (0, 1) => 1e-5, + (1, 0) => -1e-5, + (2, 3) => 2e-5, + (3, 2) => -2e-5, + (0 | 1, 2 | 3) => 1e8, + _ => 0.0, + }, + 0.0, + ) + }); + let graded = |n: usize| { + let d = |i: usize| (10.0 as Real).powf(-6.0 + 12.0 * i as Real / (n - 1) as Real); + Mat::from_fn(n, n, |i, j| sample(i, j) * scalar(d(i) / d(j), 0.0)) + }; + for a in [coupled, graded(16), graded(64), graded(128)] { + let n = a.nrows(); + let mut l = Mat::::zeros(n, n); + let mut r = Mat::::zeros(n, n); + let mut wr = Col::::zeros(n); + let mut wi = Col::::zeros(n); + let p = Spec::new(EvdParams { ..auto!(Scalar) }); + let mut mem = MemBuffer::new(evd_scratch::( + n, + ComputeEigenvectors::Yes, + ComputeEigenvectors::Yes, + Par::Seq.shape(), + p, + )); + evd_call( + a.as_ref(), + wr.as_mut(), + wi.as_mut(), + Some(l.as_mut()), + Some(r.as_mut()), + Par::Seq, + MemStack::new(&mut mem), + p, + ) + .unwrap(); + let wr: Vec<_> = wr.iter().copied().collect(); + let wi: Vec<_> = wi.iter().copied().collect(); + for (x, left) in [(l.as_ref(), true), (r.as_ref(), false)] { + let (res, _) = eigenvector_ratios(a.as_ref(), x, &wr, &wi, left); + crate::assert!( + res < get22_threshold(n), + "n {n}, left {left}: residual ratio {res}" + ); + } + } +} +/// A standardized real Schur form whose upper 2 x 2 block, shifted by the lower pair's +/// eigenvalue `i q`, has a small determinant (about `-2^-29`) next to large entries (`2^19`). +/// Raising only the determinant to a floor while keeping the block's adjugate gives a residual of +/// about `1e14`; `?LALN2` instead raises the second pivot of its elimination and solves that +/// perturbed system (#145 review). Real types only: the complex adapters take a triangular form. +#[test] +fn lapack_nonsymmetric_evd_pair_solve_with_a_small_determinant() { + if !::IS_REAL { + return; + } + let m = (2.0 as Real).powi(19); + let q = 1.0 + (2.0 as Real).powi(-30); + let t = Mat::from_fn(4, 4, |i, j| { + scalar( + match (i, j) { + (0, 1) | (0, 2) => m, + (1, 0) => -1.0 / m, + (1, 3) => -1.0, + (2, 3) => q, + (3, 2) => -q, + _ => 0.0, + }, + 0.0, + ) + }); + let mut x = Mat::::zeros(4, 4); + let p = Spec::new(EvdFromSchurParams { ..auto!(Scalar) }); + let mut mem = MemBuffer::new(eigenvectors_from_schur_scratch::(4)); + vectors_call(t.as_ref(), None, x.as_mut(), Par::Seq, MemStack::new(&mut mem), p).unwrap(); + let wr = [0.0 as Real; 4]; + let wi = [1.0, -1.0, q, -q]; + let (res, _) = eigenvector_ratios(t.as_ref(), x.as_ref(), &wr, &wi, false); + crate::assert!(res < get22_threshold(4), "residual ratio {res}"); +} +/// A 2 x 2 block `[[0, t], [-t, 0]]` above a real eigenvalue 0, with `t` below the square root of +/// the smallest positive number. Inverting the shifted block through `t^2` underflowed and gave +/// a nonfinite eigenvector (`NoConvergence`); `?LALN2`'s elimination never forms `t^2` (#145 +/// review). Real types only, as above. +#[test] +fn lapack_nonsymmetric_evd_real_shift_block_below_the_square_root_of_min_positive() { + if !::IS_REAL { + return; + } + let t = if core::mem::size_of::() == 8 { + (2.0 as Real).powi(-600) + } else { + (2.0 as Real).powi(-80) + }; + let a = Mat::from_fn(3, 3, |i, j| { + scalar( + match (i, j) { + (0, 1) => t, + (1, 0) => -t, + (0, 2) => 1.0, + _ => 0.0, + }, + 0.0, + ) + }); + let mut x = Mat::::zeros(3, 3); + let p = Spec::new(EvdFromSchurParams { ..auto!(Scalar) }); + let mut mem = MemBuffer::new(eigenvectors_from_schur_scratch::(3)); + vectors_call(a.as_ref(), None, x.as_mut(), Par::Seq, MemStack::new(&mut mem), p).unwrap(); + let (res, _) = eigenvector_ratios(a.as_ref(), x.as_ref(), &[0.0; 3], &[t, -t, 0.0], false); + crate::assert!(res < get22_threshold(3), "residual ratio {res}"); +} diff --git a/upstream-reports/README.md b/upstream-reports/README.md index e4fcd4ac..cfd26514 100644 --- a/upstream-reports/README.md +++ b/upstream-reports/README.md @@ -44,6 +44,7 @@ The program prints the expected and the actual value, and exits with status 1 wh | [faer-028](faer-028-sparse-lu-numeric-singular/) | Sparse LU factors a numerically singular matrix without an error | `faer` 0.24.4 | NaN | tensor4all/faer#132 | not yet | 2026-10-29 | | [faer-029](faer-029-evd-input-not-scaled/) | The nonsymmetric EVD loses every eigenvector digit when the entries are tiny | `faer` 0.24.4 | inaccurate value | tensor4all/faer#112 | not yet | 2026-10-29 | | [faer-030](faer-030-real-schur-aed-sort-hang/) | Real Schur QR does not return when the deflation-window sort cannot swap a pair | `faer` 0.24.4 | hang | tensor4all/faer#144 | not yet | 2026-10-29 | +| [faer-031](faer-031-real-evd-pair-clamp/) | The real nonsymmetric EVD returns a wrong eigenvector for a complex pair | `faer` 0.24.4 | wrong value | fix for tensor4all/faer#145 | not yet | 2026-10-29 | | [gemm-x86-001](gemm-x86-001-diag-negative-stride/) | A product with a diagonal gives a wrong value when the left matrix is stored backwards | `private-gemm-x86` 0.1.22 | wrong value | tensor4all/t4a-private-gemm-x86#5 (`84ca6db`) | not yet | 2026-10-29 | "Fixed upstream" was checked against upstream `main`: `faer` at `0ecd846` (2026-10-22) and `private-gemm-x86` at `71872b9` (its 0.1.22 release). diff --git a/upstream-reports/faer-031-real-evd-pair-clamp/Cargo.lock b/upstream-reports/faer-031-real-evd-pair-clamp/Cargo.lock new file mode 100644 index 00000000..5fe8ef93 --- /dev/null +++ b/upstream-reports/faer-031-real-evd-pair-clamp/Cargo.lock @@ -0,0 +1,678 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 4 + +[[package]] +name = "autocfg" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" + +[[package]] +name = "bitflags" +version = "2.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3ded4057c258ba199e2d26386d3af3780957ecaee6c4ef4041c6b4b8b97c0b06" + +[[package]] +name = "bytemuck" +version = "1.25.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" +dependencies = [ + "bytemuck_derive", +] + +[[package]] +name = "bytemuck_derive" +version = "1.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6a1f896587b6f2c069c73d2f0913e2d590c3990285cd2f0b6aa02b786b4c679c" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.9", +] + +[[package]] +name = "byteorder" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" + +[[package]] +name = "cfg-if" +version = "1.0.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4e7648175b45a9a48536d676f68d918270699102aa8dab5496df06904c914600" + +[[package]] +name = "crunchy" +version = "0.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "460fbee9c2c2f33933d720630a6a0bac33ba7053db5344fac858d4b8952d77d5" + +[[package]] +name = "defer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "930c7171c8df9fb1782bdf9b918ed9ed2d33d1d22300abb754f9085bc48bf8e8" + +[[package]] +name = "dyn-stack" +version = "0.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1c4713e43e2886ba72b8271aa66c93d722116acf7a75555cce11dcde84388fe8" +dependencies = [ + "bytemuck", + "dyn-stack-macros", +] + +[[package]] +name = "dyn-stack-macros" +version = "0.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e1d926b4d407d372f141f93bb444696142c29d32962ccbd3531117cf3aa0bfa9" + +[[package]] +name = "enum-as-inner" +version = "0.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a1e6a265c649f3f5979b601d26f1d05ada116434c87741c9493cb56218f76cbc" +dependencies = [ + "heck", + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "equator" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c35da53b5a021d2484a7cc49b2ac7f2d840f8236a286f84202369bd338d761ea" +dependencies = [ + "equator-macro 0.2.1", +] + +[[package]] +name = "equator" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "02da895aab06bbebefb6b2595f6d637b18c9ff629b4cd840965bb3164e4194b0" +dependencies = [ + "equator-macro 0.6.0", +] + +[[package]] +name = "equator-macro" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3bf679796c0322556351f287a51b49e48f7c4986e727b5dd78c972d30e2e16cc" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "equator-macro" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b14b339eb76d07f052cdbad76ca7c1310e56173a138095d3bf42a23c06ef5d8" + +[[package]] +name = "faer" +version = "0.24.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5ab6df3dd147fe8d702a288b95bcd8fcc499ab572fc80da6828f60cd4d524d67" +dependencies = [ + "bytemuck", + "dyn-stack", + "equator 0.6.0", + "faer-traits", + "gemm", + "generativity", + "libm", + "nano-gemm", + "num-complex", + "num-traits", + "private-gemm-x86", + "pulp", + "reborrow", +] + +[[package]] +name = "faer-031-real-evd-pair-clamp" +version = "0.0.0" +dependencies = [ + "faer", +] + +[[package]] +name = "faer-traits" +version = "0.24.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b87d23ed7ab1f26c0cba0e5b9e061a796fbb7dc170fa8bee6970055a1308bb0f" +dependencies = [ + "bytemuck", + "dyn-stack", + "generativity", + "libm", + "num-complex", + "num-traits", + "pulp", + "qd", + "reborrow", +] + +[[package]] +name = "gemm" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa0673db364b12263d103b68337a68fbecc541d6f6b61ba72fe438654709eacb" +dependencies = [ + "dyn-stack", + "gemm-c32", + "gemm-c64", + "gemm-common", + "gemm-f16", + "gemm-f32", + "gemm-f64", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-c32" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "086936dbdcb99e37aad81d320f98f670e53c1e55a98bee70573e83f95beb128c" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-c64" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "20c8aeeeec425959bda4d9827664029ba1501a90a0d1e6228e48bef741db3a3f" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-common" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "88027625910cc9b1085aaaa1c4bc46bb3a36aad323452b33c25b5e4e7c8e2a3e" +dependencies = [ + "bytemuck", + "dyn-stack", + "half", + "libm", + "num-complex", + "num-traits", + "once_cell", + "paste", + "pulp", + "raw-cpuid", + "seq-macro", + "sysctl", +] + +[[package]] +name = "gemm-f16" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e3df7a55202e6cd6739d82ae3399c8e0c7e1402859b30e4cb780e61525d9486e" +dependencies = [ + "dyn-stack", + "gemm-common", + "gemm-f32", + "half", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-f32" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "02e0b8c9da1fbec6e3e3ab2ce6bc259ef18eb5f6f0d3e4edf54b75f9fd41a81c" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-f64" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "056131e8f2a521bfab322f804ccd652520c79700d81209e9d9275bbdecaadc6a" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "generativity" +version = "1.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d2c81fb5260e37854d09d5c87183309fd8c555b75289427884b25660bc87a85e" + +[[package]] +name = "half" +version = "2.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6ea2d84b969582b4b1864a92dc5d27cd2b77b622a8d79306834f1be5ba20d84b" +dependencies = [ + "bytemuck", + "cfg-if", + "crunchy", + "num-traits", + "zerocopy", +] + +[[package]] +name = "heck" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" + +[[package]] +name = "hermit-abi" +version = "0.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e17592d60ebacc7d5e169f4663c5f84f9161cc90328abcfe8456f41e4dfcb284" + +[[package]] +name = "interpol" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eb58032ba748f4010d15912a1855a8a0b1ba9eaad3395b0c171c09b3b356ae50" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "libc" +version = "0.2.191" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "65cd87c9a6d60536580b316ad985726104730a5151c99845f2490428ef91240a" + +[[package]] +name = "libm" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6d2cec3eae94f9f509c767b45932f1ada8350c4bdb85af2fcab4a3c14807981" + +[[package]] +name = "nano-gemm" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9e04345dc84b498ff89fe0d38543d1f170da9e43a2c2bcee73a0f9069f72d081" +dependencies = [ + "equator 0.2.2", + "nano-gemm-c32", + "nano-gemm-c64", + "nano-gemm-codegen", + "nano-gemm-core", + "nano-gemm-f32", + "nano-gemm-f64", + "num-complex", +] + +[[package]] +name = "nano-gemm-c32" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0775b1e2520e64deee8fc78b7732e3091fb7585017c0b0f9f4b451757bbbc562" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", + "num-complex", +] + +[[package]] +name = "nano-gemm-c64" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9af49a20d58816e6b5ee65f64142e50edb5eba152678d4bb7377fcbf63f8437a" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", + "num-complex", +] + +[[package]] +name = "nano-gemm-codegen" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6cc8d495c791627779477a2cf5df60049f5b165342610eb0d76bee5ff5c5d74c" + +[[package]] +name = "nano-gemm-core" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d998dfa644de87a0f8660e5ea511d7cb5c33b5a2d9847b7af57a2565105089f0" + +[[package]] +name = "nano-gemm-f32" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "879d962e79bc8952e4ad21ca4845a21132540ed3f5e01184b2ff7f720e666523" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", +] + +[[package]] +name = "nano-gemm-f64" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9a513473dce7dc00c7e7c318481ca4494034e76997218d8dad51bd9f007a815" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", +] + +[[package]] +name = "num-complex" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" +dependencies = [ + "bytemuck", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", + "libm", +] + +[[package]] +name = "num_cpus" +version = "1.17.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91df4bbde75afed763b708b7eee1e8e7651e02d97f6d5dd763e89367e957b23b" +dependencies = [ + "hermit-abi", + "libc", +] + +[[package]] +name = "once_cell" +version = "1.21.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50" + +[[package]] +name = "paste" +version = "1.0.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" + +[[package]] +name = "private-gemm-x86" +version = "0.1.22" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bba0a7831d1297712e7e7b33759b0c837941e79196dd05df506d5f36f7f78619" +dependencies = [ + "defer", + "interpol", + "num_cpus", + "raw-cpuid", +] + +[[package]] +name = "proc-macro2" +version = "1.0.107" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "985e7ec9bb745e6ce6535b544d84d6cd6f7ad8bd711c398938ae983b91a766d9" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "pulp" +version = "0.22.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "046aa45b989642ec2e4717c8e72d677b13edd831a4d3b6cf37d9a3e54912496a" +dependencies = [ + "bytemuck", + "cfg-if", + "libm", + "num-complex", + "paste", + "pulp-wasm-simd-flag", + "raw-cpuid", + "reborrow", + "version_check", +] + +[[package]] +name = "pulp-wasm-simd-flag" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1d8f70e07b9c3962945a74e59ca1c511bba65b6419468acc217c457d93f3c740" + +[[package]] +name = "qd" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "15f1304a5aecdcfe9ee72fbba90aa37b3aa067a69d14cb7f3d9deada0be7c07c" +dependencies = [ + "bytemuck", + "libm", + "num-traits", + "pulp", +] + +[[package]] +name = "quote" +version = "1.0.47" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fbf4db142a473a8d80c26bbf18454ed458bf8d26c8219c331daecfdbd079001" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "raw-cpuid" +version = "11.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "498cd0dc59d73224351ee52a95fee0f1a617a2eae0e7d9d720cc622c73a54186" +dependencies = [ + "bitflags", +] + +[[package]] +name = "reborrow" +version = "0.5.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "03251193000f4bd3b042892be858ee50e8b3719f2b08e5833ac4353724632430" + +[[package]] +name = "same-file" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93fc1dc3aaa9bfed95e02e6eadabb4baf7e3078b0bd1b4d7b6b0b68378900502" +dependencies = [ + "winapi-util", +] + +[[package]] +name = "seq-macro" +version = "0.3.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1bc711410fbe7399f390ca1c3b60ad0f53f80e95c5eb935e52268a0e2cd49acc" + +[[package]] +name = "syn" +version = "1.0.109" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "2.0.119" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "872831b642d1a07999a962a351ed35b955ea2cfc8f3862091e2a240a84f17297" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "3.0.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d78c8dee4c7bf0e14673097256fed6142ce9d3b85a408189d07482442145823b" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "sysctl" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "01198a2debb237c62b6826ec7081082d951f46dbb64b0e8c7649a452230d1dfc" +dependencies = [ + "bitflags", + "byteorder", + "enum-as-inner", + "libc", + "thiserror", + "walkdir", +] + +[[package]] +name = "thiserror" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6aaf5339b578ea85b50e080feb250a3e8ae8cfcdff9a461c9ec2904bc923f52" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "unicode-ident" +version = "1.0.27" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a2c754d6c33795a1c324727428e5a7dedb5b06195f9890bdbcba760d3e246563" + +[[package]] +name = "version_check" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" + +[[package]] +name = "walkdir" +version = "2.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "29790946404f91d9c5d06f9874efddea1dc06c5efe94541a7d6863108e3a5e4b" +dependencies = [ + "same-file", + "winapi-util", +] + +[[package]] +name = "winapi-util" +version = "0.1.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c2a7b1c03c876122aa43f3020e6c3c3ee5c05081c9a00739faf7503aeba10d22" +dependencies = [ + "windows-sys", +] + +[[package]] +name = "windows-link" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f0805222e57f7521d6a62e36fa9163bc891acd422f971defe97d64e70d0a4fe5" + +[[package]] +name = "windows-sys" +version = "0.61.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ae137229bcbd6cdf0f7b80a31df61766145077ddf49416a728b02cb3921ff3fc" +dependencies = [ + "windows-link", +] + +[[package]] +name = "zerocopy" +version = "0.8.63" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e5fe1f8f1b06191a00962174c61aa5005e0bb391a6d80d07e24d115c01a92ed8" +dependencies = [ + "zerocopy-derive", +] + +[[package]] +name = "zerocopy-derive" +version = "0.8.63" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "863ad3ac83293fb4d740aedbfdc9240dd8d1a50c1099acd76ce80ce7c7230c7f" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] diff --git a/upstream-reports/faer-031-real-evd-pair-clamp/Cargo.toml b/upstream-reports/faer-031-real-evd-pair-clamp/Cargo.toml new file mode 100644 index 00000000..9aad4a0b --- /dev/null +++ b/upstream-reports/faer-031-real-evd-pair-clamp/Cargo.toml @@ -0,0 +1,11 @@ +[package] +name = "faer-031-real-evd-pair-clamp" +version = "0.0.0" +edition = "2021" +publish = false + +# Not part of the t4a-faer workspace: this reproduces the published crate on its own. +[workspace] + +[dependencies] +faer = { version = "=0.24.4", default-features = false, features = ["std", "linalg"] } diff --git a/upstream-reports/faer-031-real-evd-pair-clamp/README.md b/upstream-reports/faer-031-real-evd-pair-clamp/README.md new file mode 100644 index 00000000..8d613027 --- /dev/null +++ b/upstream-reports/faer-031-real-evd-pair-clamp/README.md @@ -0,0 +1,56 @@ +# faer-031: the real nonsymmetric EVD returns a wrong eigenvector for a complex pair + +**What happens:** `Mat::eigen` of a real matrix returns an eigenvector of a complex eigenvalue pair whose residual `|A u - lambda u|` is about 9e4 `eps` of `|A| |u|`, where LAPACK's `dgeev` gives 1e-13. +**Which input:** a real matrix with two complex pairs that are small and close compared with the rest of the matrix. The program below uses two 2 x 2 rotation blocks (eigenvalues `+-1e-5 i` and `+-2e-5 i`) coupled by entries of size `1e8`. +**Which version:** `faer` 0.24.4. +**Added to this list:** 2026-10-29. + +## Reproduce + +```sh +cargo run +``` + +The program builds the 4 x 4 matrix shown at the top of `src/main.rs`. +It calls `Mat::eigen` and prints `|A u - lambda u|_1 / (|A|_1 |u|_1 eps)` for each eigenvector. + +| eigenvalue | expected (LAPACK `dgeev`) | `faer` 0.24.4 | +|---|---|---| +| `+-1e-5 i` | about 1e-13 | 0 | +| `+-2e-5 i` | about 1e-13 | **8.97e4** | + +## Cause + +The eigenvectors of a real matrix come from a back substitution on its real Schur form `T`. +For a complex eigenvalue `lambda = p + i q`, each 2 x 2 diagonal block `B` of `T` above it is solved as `(B - lambda I) x = r` (a 2 x 2 complex system). +`solve_cplx_shifted_upper_quasi_triangular_system` in `src/linalg/evd/mod.rs` solves it with the inverse of the shifted block, through its determinant `z = det(B - lambda I)` (line 825). +It guards against a nearly singular block with `if z < eps * norm` (line 828). +Here `norm` is the sum of the magnitudes of all entries of `T` (lines 567 to 570), so `eps * norm` has the size of an entry. +A determinant is a product of two entries, so it is compared with a threshold of the wrong size. +In this matrix the shifted blocks have entries near `1e-5` and determinants near `1e-10`, below `eps * norm` (about `9e-8`, from the entries `1e8`). +The determinant is then replaced by `eps * norm` while the numerators of the inverse stay those of the original block, so the two entries of that row are solved from a system that does not exist, and they are wrong. +The 1 x 1 case in the same function (line 794) compares `|T[i, i] - lambda|` with `eps * norm`, which has the right size. + +LAPACK's `dtrevc3` solves the same blocks with `dlaln2`. +That routine compares its threshold with the largest entry of the shifted block. +It then solves the block by Gaussian elimination with complete pivoting, and raises the second pivot to the threshold when it is smaller, so it always solves a slightly perturbed system exactly. + +## How the fork fixed it + +The fork solves each shifted 2 x 2 block the way `dlaln2` does. +A block whose entries are all below the threshold is taken as the threshold times the identity. +Any other block is solved by elimination with complete pivoting, with the second pivot raised to the threshold when it is smaller. +The 2 x 2 blocks above a real eigenvalue use the same solve with a zero imaginary part. +With that solve in place of lines 825 to 843 of the 0.24.4 source, this program prints residual ratios of 7e-14 for all four eigenvectors. + +Fork regression tests: `linalg::evd::lapack::{f32,f64}::lapack_nonsymmetric_evd_graded_pairs_keep_their_residual` (this matrix, and strongly graded matrices `D A D^-1` with `D = diag(10^linspace(-6, 6, n))`, left and right eigenvectors) and `lapack_nonsymmetric_evd_pair_solve_with_a_small_determinant`. + +## Affected versions + +Confirmed by running this program: `faer` 0.24.4. +Not fixed upstream: the same code is in upstream `main` at `0ecd846` (2026-10-22). + +## References + +* Fork fix: the pull request that closes tensor4all/faer#145. +* Where it was found: tensor4all/faer#145, strongly graded matrices in the fork's comparison with LAPACK. In the fork, a Schur form is scaled to entries near one before this solve, which made every 2 x 2 block of such a matrix hit the guard. diff --git a/upstream-reports/faer-031-real-evd-pair-clamp/src/main.rs b/upstream-reports/faer-031-real-evd-pair-clamp/src/main.rs new file mode 100644 index 00000000..9d55dfb4 --- /dev/null +++ b/upstream-reports/faer-031-real-evd-pair-clamp/src/main.rs @@ -0,0 +1,39 @@ +// A 4 x 4 real matrix with two complex eigenvalue pairs, +-1e-5 i and +-2e-5 i, +// in two 2 x 2 rotation blocks that are coupled by entries of size 1e8: +// +// [ 0 1e-5 1e8 1e8 ] +// [ -1e-5 0 1e8 1e8 ] +// [ 0 0 0 2e-5 ] +// [ 0 0 -2e-5 0 ] +// +// We check every eigenvector u with eigenvalue lambda through the residual ratio +// |A u - lambda u|_1 / (|A|_1 |u|_1 eps). LAPACK's dgeev gives 9.6e-14 here. +use faer::prelude::*; + +fn main() { + let a = Mat::::from_fn(4, 4, |i, j| match (i, j) { + (0, 1) => 1e-5, + (1, 0) => -1e-5, + (2, 3) => 2e-5, + (3, 2) => -2e-5, + (0 | 1, 2 | 3) => 1e8, + _ => 0.0, + }); + let evd = a.eigen().unwrap(); + let (u, s) = (evd.U(), evd.S()); + + let ac = Mat::::from_fn(4, 4, |i, j| c64::new(a[(i, j)], 0.0)); + let col_norm = |m: MatRef<'_, c64>, j: usize| (0..4).map(|i| m[(i, j)].norm()).sum::(); + let a_norm = (0..4).map(|j| col_norm(ac.as_ref(), j)).fold(0.0, f64::max); + let r = &ac * u - u * s.as_ref(); + let mut worst = 0.0f64; + for j in 0..4 { + let ratio = col_norm(r.as_ref(), j) / (a_norm * col_norm(u, j) * f64::EPSILON); + println!("eigenvalue {:+.1e}i: residual ratio {ratio:.2e}", s[j].im); + worst = worst.max(ratio); + } + if worst > 20.0 { + println!("WRONG: an eigenvector is far from exact (dgeev: 9.6e-14)"); + std::process::exit(1); + } +} From 3ea3887f6646e347d40f35924a5c6735fd00d804 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Thu, 29 Oct 2026 21:04:48 +0900 Subject: [PATCH 6/6] PATCHSET: the nonsymmetric EVD fixes land together (#112, #144, #145) With #112's input scaling in the same change the c32 known failure is gone, so the #145 row no longer reports it; it records the combined result. Co-Authored-By: Claude Opus 5.5 --- PATCHSET.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PATCHSET.md b/PATCHSET.md index fedee0bc..ea5ca014 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -102,7 +102,7 @@ gate has changed and two kinds of adoption are now in the ledger. | 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/). | -| The 2 x 2 block solves of the real eigenvector back substitution follow `?LALN2` (#145) | `solve_cplx_shifted_upper_quasi_triangular_system` solved each 2 x 2 diagonal block of the real Schur form for a complex eigenvalue through the block's inverse, and compared its determinant with the entry-sized clamp `eps_norm`. The determinant has the size of a product of two entries, so every block with entries below `sqrt(eps_norm)` got the determinant `eps_norm` with its original numerators. After `scaled_into_range` (#112) the blocks of a strongly graded matrix `D A D^-1` (`D = diag(10^linspace(-6, 6, n))`) are all that small, and every complex pair got a wrong eigenvector: GET22 `1.0e9` at `f64` n = 64 in the issue's sweep. Stage comparison against LAPACK 3.12.0 (`f64` n = 64): the reduction is fine (Schur backward error 8.0 eps |A|); on faer's own Schur form `dtrevc3` gives `2.2e-4` and faer `3.4e7`, and all 24 failing eigen-blocks are complex pairs. Both 2 x 2 solves (complex shift, and the real-shift one above a real eigenvalue, which inverted through `a^2 - b c` and gave `inf` when that underflowed) now use `solve_shifted_2x2_cplx`, `?LALN2`'s method: `eps_norm * I` when every entry is below `eps_norm`, otherwise elimination with complete pivoting and the second pivot raised to `eps_norm`. Raising only the determinant was tried first and rejected in review: with the original numerators it is not a solved system (residual `4.5e14` on a 4 x 4 Schur form with a small determinant next to large entries). Balancing (`?GEBAL`) is not the cause: `dgeev` without balancing gives 0.58 on the same matrix, the fixed faer 0.45; balancing would bring both to about `1e-10`. Upstream has the same comparison (`z < eps * norm`, faer 0.24.4 `src/linalg/evd/mod.rs:828`); without the fork's scaling it fails on smaller inputs, for example two small, close complex pairs coupled by large entries. | KEEP | tensor4all/faer#145 | Correctness, no timing claim: one 2 x 2 elimination per block in place of one inverse. Tests (four dtypes, the real ones reach the changed code): `lapack_nonsymmetric_evd_graded_pairs_keep_their_residual` (graded n = 16, 64, 128 and the coupled 4 x 4, left and right, GET22 against `get22_threshold`; before: `5.9e8` at n = 16, `1.5e12` on the 4 x 4), `lapack_nonsymmetric_evd_pair_solve_with_a_small_determinant` (fails with the determinant floor: `4.5e14`; on `main`: `f32` `4.2e6`) and `lapack_nonsymmetric_evd_real_shift_block_below_the_square_root_of_min_positive` (`NoConvergence` on `main`). `linalg::evd` with `--include-ignored`: the one failure is the known c32 kinds 6/7 case of #112. Upstream report: [`upstream-reports/faer-031-real-evd-pair-clamp`](upstream-reports/faer-031-real-evd-pair-clamp/). | +| The 2 x 2 block solves of the real eigenvector back substitution follow `?LALN2` (#145) | `solve_cplx_shifted_upper_quasi_triangular_system` solved each 2 x 2 diagonal block of the real Schur form for a complex eigenvalue through the block's inverse, and compared its determinant with the entry-sized clamp `eps_norm`. The determinant has the size of a product of two entries, so every block with entries below `sqrt(eps_norm)` got the determinant `eps_norm` with its original numerators. After `scaled_into_range` (#112) the blocks of a strongly graded matrix `D A D^-1` (`D = diag(10^linspace(-6, 6, n))`) are all that small, and every complex pair got a wrong eigenvector: GET22 `1.0e9` at `f64` n = 64 in the issue's sweep. Stage comparison against LAPACK 3.12.0 (`f64` n = 64): the reduction is fine (Schur backward error 8.0 eps |A|); on faer's own Schur form `dtrevc3` gives `2.2e-4` and faer `3.4e7`, and all 24 failing eigen-blocks are complex pairs. Both 2 x 2 solves (complex shift, and the real-shift one above a real eigenvalue, which inverted through `a^2 - b c` and gave `inf` when that underflowed) now use `solve_shifted_2x2_cplx`, `?LALN2`'s method: `eps_norm * I` when every entry is below `eps_norm`, otherwise elimination with complete pivoting and the second pivot raised to `eps_norm`. Raising only the determinant was tried first and rejected in review: with the original numerators it is not a solved system (residual `4.5e14` on a 4 x 4 Schur form with a small determinant next to large entries). Balancing (`?GEBAL`) is not the cause: `dgeev` without balancing gives 0.58 on the same matrix, the fixed faer 0.45; balancing would bring both to about `1e-10`. Upstream has the same comparison (`z < eps * norm`, faer 0.24.4 `src/linalg/evd/mod.rs:828`); without the fork's scaling it fails on smaller inputs, for example two small, close complex pairs coupled by large entries. | KEEP | tensor4all/faer#145 | Correctness, no timing claim: one 2 x 2 elimination per block in place of one inverse. Tests (four dtypes, the real ones reach the changed code): `lapack_nonsymmetric_evd_graded_pairs_keep_their_residual` (graded n = 16, 64, 128 and the coupled 4 x 4, left and right, GET22 against `get22_threshold`; before: `5.9e8` at n = 16, `1.5e12` on the 4 x 4), `lapack_nonsymmetric_evd_pair_solve_with_a_small_determinant` (fails with the determinant floor: `4.5e14`; on `main`: `f32` `4.2e6`) and `lapack_nonsymmetric_evd_real_shift_block_below_the_square_root_of_min_positive` (`NoConvergence` on `main`). Merged in one change with #112's input scaling and #144's sort: every LAPACK validation family passes (nonsymmetric EVD 40/40), and on the #112 sweep the real graded matrices (f32/f64, n = 16 to 257) give GET22 ratios of 0.20 to 0.87, down from up to 2.6e9. Upstream report: [`upstream-reports/faer-031-real-evd-pair-clamp`](upstream-reports/faer-031-real-evd-pair-clamp/). | | 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). |