diff --git a/PATCHSET.md b/PATCHSET.md index b2cb793c..584899f3 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 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 caa6e546..0a235151 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: 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 7aba8438..8347365f 100644 --- a/faer/src/linalg/evd/mod.rs +++ b/faer/src/linalg/evd/mod.rs @@ -1833,6 +1833,47 @@ pub fn eigenvectors_from_schur_cplx( Err(EvdError::NoConvergence) } } +/// `?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 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)].real().abs()).fmax(A[(i, j)].imag().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 +1920,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 +2059,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 +3154,134 @@ 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. 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; + + /// 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(); + // `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)); + } + } + } + }}; + } + + /// 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); + 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": [] } 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..b9f93f1d --- /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, 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. + +## Reproduce + +```sh +cargo run +``` + +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 | **`U` not finite, returned with `Ok`** | +| f32 | 1 | 15.8 | 10.4 | +| f32 | `2^-100` | 15.8 | **`U` not finite, returned with `Ok`** | +| f32 | `2^-110` | 15.8 | **`U` not finite, returned with `Ok`** | + +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 + +`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 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. + +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..359420e4 --- /dev/null +++ b/upstream-reports/faer-029-evd-input-not-scaled/src/main.rs @@ -0,0 +1,52 @@ +//! 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) + }; + // `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); + } + if failed { + println!("the eigenvectors of a matrix with tiny entries have lost every digit"); + std::process::exit(1); + } +}