diff --git a/PATCHSET.md b/PATCHSET.md index b2cb793c..bb13d9f5 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 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-029-real-evd-pair-clamp`](upstream-reports/faer-029-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 7aba8438..5cab5353 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 75eed75a..941f4fc8 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-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-029-real-evd-pair-clamp/Cargo.lock b/upstream-reports/faer-029-real-evd-pair-clamp/Cargo.lock new file mode 100644 index 00000000..e19f9e54 --- /dev/null +++ b/upstream-reports/faer-029-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-029-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-029-real-evd-pair-clamp/Cargo.toml b/upstream-reports/faer-029-real-evd-pair-clamp/Cargo.toml new file mode 100644 index 00000000..585ec3b8 --- /dev/null +++ b/upstream-reports/faer-029-real-evd-pair-clamp/Cargo.toml @@ -0,0 +1,11 @@ +[package] +name = "faer-029-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-029-real-evd-pair-clamp/README.md b/upstream-reports/faer-029-real-evd-pair-clamp/README.md new file mode 100644 index 00000000..63a050ee --- /dev/null +++ b/upstream-reports/faer-029-real-evd-pair-clamp/README.md @@ -0,0 +1,56 @@ +# faer-029: 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-029-real-evd-pair-clamp/src/main.rs b/upstream-reports/faer-029-real-evd-pair-clamp/src/main.rs new file mode 100644 index 00000000..9d55dfb4 --- /dev/null +++ b/upstream-reports/faer-029-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); + } +}