diff --git a/PATCHSET.md b/PATCHSET.md index 899f0fd3..b2cb793c 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -107,6 +107,7 @@ gate has changed and two kinds of adoption are now in the ledger. | Full pivoting LU: the pivot search's ties are independent of the SIMD lane count (#102, PR #105) | `reduce_2d`, the SIMD pivot search of `lu::full_pivoting`, took the first lane that held the maximum (`first_true_mask`). Which lane comes first depends on the lane count, so for a matrix with tied maxima (small integers tie a lot) an AVX-512 build and an AVX2 build chose different pivots and eliminated in a different order. On a singular fixture one of them ended on an exact zero pivot, which the solve turned into `NaN`. The tie is now broken as the scalar search breaks it: the smallest column, then the smallest row. Upstream has the same code: faer 0.24.4 `src/linalg/lu/full_pivoting/factor.rs:64` (`simd.first_true_mask(is_best)`). Not reported upstream. | KEEP | tensor4all/faer#102 | Correctness: the pivots are the same on every vector width. Test `pivot_search_breaks_ties_by_column_then_row` (the SIMD and scalar searches agree on tied matrices of several shapes). Upstream report: [`upstream-reports/faer-001-full-piv-lu-lane-tie`](upstream-reports/faer-001-full-piv-lu-lane-tie/). | | Partial SVD: a retained Ritz pair carries its residual estimate through the restart (#102, PR #108) | In the restart of the partial SVD (`operator::svd`), a retained pair that moves to a lower slot moved its value and its vector but not its residual estimate. The coupling row of the next iteration is built from the residual array (`H[(i, k)] = residual[i]`), so the moved pair coupled with another pair's estimate. The locked branch already moved it. The retained branch now moves it too. Upstream has the same code: faer 0.24.4 `src/operator/svd/mod.rs:273-275`, whose `Group::Retain` branch moves `H[(mi, mi)]` and the vector, but not `residual`. The partial self-adjoint eigensolver (`operator::self_adjoint_eigen`) is not affected, because it rebuilds the coupling row from `H Q`. Not reported upstream. | KEEP | tensor4all/faer#102 | Correctness. PR #108 adds no dedicated test; its record is the PR. Upstream report: [`upstream-reports/faer-002-partial-svd-residual`](upstream-reports/faer-002-partial-svd-residual/). | | Krylov solvers: final residual and BiCGSTAB breakdown (#133, #134; upstream #328, #327) | When `bicgstab` or `conjugate_gradient` stopped at `max_iters`, `NoConvergence` carried the initial residual `|b|`: `bicgstab` kept each new norm in a local, and `conjugate_gradient` shadowed `abs_residual` inside the loop. Both now report the residual of the returned `x`. `bicgstab`'s rank test on `r̃ᵀ v` took its tolerance from its own largest pivot, so a zero first pivot passed and the step divided by zero, turning `x` into NaN until `max_iters`. The tolerance is now `eps * k * max_j |r̃_j| * max_j |v_j|`, the bound on one entry of `r̃ᵀ v`, so a wide block of right-hand sides is not scaled by its width. The comparison is strict, so `v = 0` also counts. At a breakdown (rank 0), or when the half step `s` or the new residual is not finite, the iteration stops. It keeps the last iterate whose residual is finite and reports that residual as `NoConvergence`. A finite residual does not bound `x` itself: an operator scaled near the range limit can still overflow `x`. `InitialGuessStatus::Zero` is now written out as zeros on entry in both solvers, so every exit, including `max_iters = 0`, returns zeros or an iterate, never the caller's buffer. A breakdown-specific error variant is the feature in #135 (upstream #334) and is not part of this change. Written from the upstream reports (Codeberg sarah-quinones/faer#327, #328); no code was taken from the open upstream PRs #330 and #331. | KEEP | tensor4all/faer#133, tensor4all/faer#134 | Correctness only, no timing claim: the rank test adds `2k` column norms of length `n` per iteration, and a zero guess adds one fill. Nine focused tests: `codeberg_328` for both solvers, `codeberg_327`, a first-step breakdown with a NaN-filled zero guess, `A r̃ = 0`, a non-finite second half step, a zero guess with `max_iters = 0` in both solvers, and an `f32` identity with 64 right-hand sides. Before the fix, `codeberg_327`, both `codeberg_328` tests and the first-step breakdown failed. All 23 `matrix_free` tests pass. Independent review: gpt-6.1-sol. | +| Sparse LU reports an exact zero pivot (#132; upstream #339, #332) | The numeric sparse LU returned `Ok` for a numerically singular matrix, and the solve returned NaN. The simplicial backend's pivot search starts from `-1`, so an exact-zero candidate wins and is inverted; `SymbolicSingular` was raised only when a column had no candidate. The supernodal backend factors each panel with the dense partial-pivoting kernel, which gives no singular signal. Upstream reports the simplicial backend only; the supernodal one has the same defect. Both backends now return the new `LuError::NumericSingular { index }` when a pivot is exactly zero: the simplicial one on the chosen pivot, the supernodal one on a zero diagonal entry of the factored panel. That is the variant the upstream maintainer asked for in upstream PR #332, closed unmerged. `index` counts elimination steps in the fill-reducing column order, as `SymbolicSingular`'s does. The new variant is an API change on an exhaustive public enum; no tensor4all consumer references the sparse `LuError`. Tolerance-based detection is #146. Written from the upstream report; no code was taken from upstream PR #332. Upstream report: [`upstream-reports/faer-028-sparse-lu-numeric-singular`](upstream-reports/faer-028-sparse-lu-numeric-singular/) (reproduced on `faer` 0.24.4, both backends). | KEEP | tensor4all/faer#132, tensor4all/faer#146 | Correctness only, no timing claim: one comparison per simplicial column and one per supernodal panel diagonal entry. Three focused tests, each in both backends: `codeberg_339` (the upstream reproduction at `f32`/`f64`/`c32`/`c64`), a 4x4 singular/nonsingular pair, and a nonzero subnormal complex pivot. The simplicial check reads the pivot itself (`abs1`), not its modulus, because the complex modulus of a nonzero subnormal underflows to zero. The first two tests fail without the checks, and all three pass with them. All 23 `sparse::linalg` tests pass. Independent review: gpt-6.1-sol. | | Tridiagonal EVD at the range limits | Both tridiagonal backends now rescale `(diag, offdiag)` when the largest entry leaves the window in which an entry squared and a reciprocal both stay finite (`[sqrt(min_positive)/eps, sqrt(max_positive)]`), and scale the computed eigenvalues back. The QR fallback used to multiply by `1/anorm`, which is `inf` for a subnormal spectrum; divide-and-conquer did not rescale at all and returned `inf`/`NaN` for a spectrum near the top of the range. The 2x2 shortcut, the eigenvectors, and the public results for an ordinary spectrum are unchanged. | KEEP | tensor4all/faer#57 | Focused tests `linalg::evd::tridiagonal_range_tests::*`: the two reported matrices, and both backends against a closed form at both ends of the range. Record: `docs/worklogs/2026-10-25-tridiagonal-evd-range.md`. Upstream report: [`upstream-reports/faer-017-tridiag-evd-range`](upstream-reports/faer-017-tridiag-evd-range/). | | Schur eigenvector extraction at the range limit | (The clamp described here is superseded by the #112 row above: it is now relative to the eigenvalue and the largest entry, on a Schur form scaled into range.) The shift clamp of the shifted back substitution is now `eps_norm = eps * norm`, computed once, where `norm` is the sum of the entry magnitudes; it is recomputed term by term (`eps * abs(re) + eps * abs(im)` per entry) only when the plain sum is not finite. The plain sum overflows to `inf` for a finite Schur form near the top of the range, which clamps every shift to `inf` and returns a wrong direction with `Ok`. Scaling the terms unconditionally was rejected because it rounds the clamp of a subnormal matrix to zero. The 2x2 block solves keep their arithmetic: at the top of the range they fail with `NoConvergence` instead of returning a wrong vector (a power-of-two rescaling loop was tried and dropped after review: it could hang on an overflowed shifted diagonal and lost the small entry of a wide-range block). `q` needs no change (geometric mean of two finite entries). | KEEP | tensor4all/faer#77 | Correctness only. Focused tests `linalg::evd::range_limit_vector_tests::*` at `f32`/`f64`/`c32`/`c64`, at scale 1 and at `2^127`/`2^1023`: the reported 2x2 fails at base with the reported 2/3 residual; conjugate-pair blocks at the range limit must fail or be correct, never a wrong `Ok`; a wide-range pair keeps its exact solution; the subnormal control fails under the unconditional term-by-term variant. The ignored LAPACK nonsymmetric/Schur run is byte-identical to base (same `LAPACK_CASES`, 1286 identical diagnostic lines, same failing set), and the library suite passes 464. Limits left: the complex reciprocal in `faer-traits` still returns zero for a shift near `2^1023`, and a range-safe 2x2 block solve is not done. Record: `docs/worklogs/2026-10-28-schur-eigenvector-range.md`. Upstream report: [`upstream-reports/faer-025-evd-norm-sum-overflow`](upstream-reports/faer-025-evd-norm-sum-overflow/). | | Pivoted LLT zero first pivot | `llt_pivoting::factor::cholesky_in_place` compared only pivots after the first with its rank tolerance `n * eps * max(diag)`, which is zero for a matrix whose diagonal is zero. The 1x1 zero matrix then returned rank one, and a larger zero matrix took a zero square root, scaled by its infinite reciprocal and returned `NonPositivePivot`. A zero largest diagonal entry now returns rank zero with a zero stored factor, as LAPACK `dpstf2`/`dpstrf` do (`AJJ.LE.ZERO` → `RANK = 0` before the pivot loop). Positive and negative diagonals are unchanged. The `j > 0` guard is upstream faer code (`64280357`); nothing was reported upstream. | KEEP | tensor4all/faer#31 | Correctness only. Found by the #31 LAPACK validation harness (`lapack_cholesky_cases`: 26 diagnostics per type, 0 after the change at all four types). Test `test_zero_matrix_rank_zero` (n = 1, 2, 5; rank and reconstruction) fails before the change and passes after; library suite 247 pass, with the 3 pre-existing strict-triangle matmul failures that `main` also has. Upstream report: [`upstream-reports/faer-019-pivoted-llt-zero-matrix`](upstream-reports/faer-019-pivoted-llt-zero-matrix/). | diff --git a/faer/src/sparse/linalg/lu.rs b/faer/src/sparse/linalg/lu.rs index 2e9eab76..bad16ce7 100644 --- a/faer/src/sparse/linalg/lu.rs +++ b/faer/src/sparse/linalg/lu.rs @@ -945,6 +945,13 @@ pub mod supernodal { par.reborrow(), params, ); + for i in 0..s_size { + if s_L[(i, i)].abs1() == zero::() { + return Err(LuError::NumericSingular { + index: s_begin + i, + }); + } + } for (idx, t) in transpositions.iter().enumerate() { let i_t = s_row_idx[idx + t.zx()].zx(); let kk = row_perm_inv[i_t].zx(); @@ -1749,6 +1756,11 @@ pub mod simplicial { return Err(LuError::SymbolicSingular { index: j }); } let x_piv = x[pivot_idx].copy(); + // the pivot itself, not `pivot_val`: a complex modulus can + // underflow to zero for a nonzero subnormal entry + if x_piv.abs1() == zero::() { + return Err(LuError::NumericSingular { index: j }); + } let ref x_piv_inv = x_piv.recip(); row_perm_inv[pivot_idx] = I(j); lu.u_row_idx[u_pos] = I(j); @@ -2765,4 +2777,119 @@ mod tests { } } } + fn factorize_lu_error( + A: SparseColMatRef<'_, usize, T>, + supernodal_flop_ratio_threshold: SupernodalThreshold, + ) -> Result<(), LuError> { + let symbolic = factorize_symbolic_lu( + A.symbolic(), + LuSymbolicParams { + supernodal_flop_ratio_threshold, + ..Default::default() + }, + ) + .unwrap(); + let mut numeric = NumericLu::::new(); + symbolic + .factorize_numeric_lu( + &mut numeric, + A, + Par::Seq, + MemStack::new(&mut MemBuffer::new( + symbolic.factorize_numeric_lu_scratch::( + ParShape::Seq, + Default::default(), + ), + )), + Default::default(), + ) + .map(|_| ()) + } + fn numeric_singular_index( + n: usize, + entries: &[(usize, usize, f64)], + supernodal_flop_ratio_threshold: SupernodalThreshold, + ) -> Option { + let triplets = entries + .iter() + .map(|&(i, j, v)| Triplet::new(i, j, from_f64::(v))) + .collect::>(); + let A = + SparseColMat::::try_new_from_triplets(n, n, &triplets) + .unwrap(); + match factorize_lu_error(A.as_ref(), supernodal_flop_ratio_threshold) { + Err(LuError::NumericSingular { index }) => Some(index), + Ok(()) => None, + Err(e) => panic!("unexpected error {e:?}"), + } + } + #[test] + fn codeberg_339() { + // upstream sarah-quinones/faer#339: `[[1, 1], [1, 1]]` factored as `Ok` + // and the solve returned NaN; both backends must report the zero pivot, + // which only the second step can meet + let A = [(0, 0, 1.0), (0, 1, 1.0), (1, 0, 1.0), (1, 1, 1.0)]; + for threshold in [ + SupernodalThreshold::FORCE_SIMPLICIAL, + SupernodalThreshold::FORCE_SUPERNODAL, + ] { + assert!(numeric_singular_index::(2, &A, threshold) == Some(1)); + assert!(numeric_singular_index::(2, &A, threshold) == Some(1)); + assert!(numeric_singular_index::(2, &A, threshold) == Some(1)); + assert!(numeric_singular_index::(2, &A, threshold) == Some(1)); + } + } + #[test] + fn numeric_singular_after_an_elimination_step() { + // the second column is a multiple of the first, so one elimination step + // meets a zero pivot. `index` counts steps in the fill-reducing column + // order, as `SymbolicSingular`'s does, so only its range is checked + let A = [ + (0, 0, 2.0), + (1, 0, 4.0), + (0, 1, 1.0), + (1, 1, 2.0), + (2, 2, 1.0), + (0, 3, 1.0), + (3, 3, 1.0), + (2, 3, 1.0), + ]; + // a structurally full matrix that is nonsingular still factors + let B = [ + (0, 0, 2.0), + (1, 0, 4.0), + (0, 1, 1.0), + (1, 1, 3.0), + (2, 2, 1.0), + (0, 3, 1.0), + (3, 3, 1.0), + (2, 3, 1.0), + ]; + for threshold in [ + SupernodalThreshold::FORCE_SIMPLICIAL, + SupernodalThreshold::FORCE_SUPERNODAL, + ] { + assert!(matches!( + numeric_singular_index::(4, &A, threshold), + Some(index) if index < 4 + )); + assert!(numeric_singular_index::(4, &B, threshold) == None); + } + } + #[test] + fn nonzero_subnormal_pivot_is_not_singular() { + // the modulus of `c64(2^-1074, 0)` underflows to zero, the entry + // itself is not zero + let tiny = f64::from_bits(1); + let triplets = [Triplet::new(0, 0, c64::new(tiny, 0.0))]; + let A = + SparseColMat::::try_new_from_triplets(1, 1, &triplets) + .unwrap(); + for threshold in [ + SupernodalThreshold::FORCE_SIMPLICIAL, + SupernodalThreshold::FORCE_SUPERNODAL, + ] { + assert!(factorize_lu_error(A.as_ref(), threshold).is_ok()); + } + } } diff --git a/faer/src/sparse/linalg/mod.rs b/faer/src/sparse/linalg/mod.rs index c57f411c..006aa381 100644 --- a/faer/src/sparse/linalg/mod.rs +++ b/faer/src/sparse/linalg/mod.rs @@ -77,6 +77,12 @@ pub enum LuError { /// iteration at which a pivot could not be found index: usize, }, + /// rank deficient numeric values: every candidate pivot of a column is + /// exactly zero + NumericSingular { + /// iteration at which the pivot is zero + index: usize, + }, /// non algorithmic error Generic(FaerError), } diff --git a/upstream-reports/README.md b/upstream-reports/README.md index 0d71e03f..75eed75a 100644 --- a/upstream-reports/README.md +++ b/upstream-reports/README.md @@ -4,6 +4,7 @@ This list is updated over time. Last updated: 2026-10-29. This directory lists defects that we found while working on a port of `faer` and that also occur in the published upstream crates. Each entry was reproduced by running its program against the upstream release named in the table. +faer-028 was first reported upstream (sarah-quinones/faer#339, simplicial backend only). It is listed because the supernodal backend has the same defect. Defects that exist only in this fork are not listed here. Each subdirectory is a standalone Cargo project. @@ -40,6 +41,7 @@ The program prints the expected and the actual value, and exits with status 1 wh | [faer-025](faer-025-evd-norm-sum-overflow/) | Nonsymmetric EVD returns a wrong eigenvector when the Schur norm sum overflows | `faer` 0.24.4 | wrong value | tensor4all/faer#121 (`0c806f1`) | not yet | 2026-10-29 | | [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 | | [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-028-sparse-lu-numeric-singular/Cargo.lock b/upstream-reports/faer-028-sparse-lu-numeric-singular/Cargo.lock new file mode 100644 index 00000000..ab6a1e38 --- /dev/null +++ b/upstream-reports/faer-028-sparse-lu-numeric-singular/Cargo.lock @@ -0,0 +1,679 @@ +# 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-028-sparse-lu-numeric-singular" +version = "0.0.0" +dependencies = [ + "dyn-stack", + "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-028-sparse-lu-numeric-singular/Cargo.toml b/upstream-reports/faer-028-sparse-lu-numeric-singular/Cargo.toml new file mode 100644 index 00000000..d0877dea --- /dev/null +++ b/upstream-reports/faer-028-sparse-lu-numeric-singular/Cargo.toml @@ -0,0 +1,12 @@ +[package] +name = "faer-028-sparse-lu-numeric-singular" +version = "0.0.0" +edition = "2021" +publish = false + +# A standalone project: it reproduces the published crate on its own. +[workspace] + +[dependencies] +faer = { version = "=0.24.4", default-features = false, features = ["std", "linalg", "sparse"] } +dyn-stack = "0.13" diff --git a/upstream-reports/faer-028-sparse-lu-numeric-singular/README.md b/upstream-reports/faer-028-sparse-lu-numeric-singular/README.md new file mode 100644 index 00000000..704fee24 --- /dev/null +++ b/upstream-reports/faer-028-sparse-lu-numeric-singular/README.md @@ -0,0 +1,46 @@ +# faer-028: sparse LU factors a numerically singular matrix without an error + +**What happens:** the numeric sparse LU factorization returns `Ok` for a matrix whose pivot is exactly zero, and the solve then returns NaN. +**Which input:** a structurally full but singular matrix, in both the simplicial and the supernodal backend. The program uses `A = [[1, 1], [1, 1]]`. +**Which version:** `faer` 0.24.4. +**Added to this list:** 2026-10-29. + +## Reproduce + +```sh +cargo run +``` + +The program forces each backend in turn (`SupernodalThreshold::FORCE_SIMPLICIAL`, then `FORCE_SUPERNODAL`), factors `A`, and, if the factorization succeeds, solves `A x = [1, 1]`. + +| backend | expected | `faer` 0.24.4 | +|---|---|---| +| simplicial | an error | **`Ok`, then `x = [NaN, NaN]`** | +| supernodal | an error | **`Ok`, then `x = [NaN, NaN]`** | + +## Cause + +- **Simplicial.** The pivot search starts from `pivot_val = -1` (`src/sparse/linalg/lu.rs:1729`), so any candidate wins, including an exact zero. `SymbolicSingular` is raised only when a column has no candidate at all (line 1748). Nothing rejects the zero pivot: it is stored in `U`, and its reciprocal (line 1752) scales whatever rows remain below it. +- **Supernodal.** Each supernode's panel is factored by the dense `lu_in_place_recursion` (line 940), which gives no singular signal. + +In this program the zero pivot is the last one, so no rows remain below it, and the factors stay finite with a zero on `U`'s diagonal. The solve divides by that zero, which gives the NaN. When rows do remain below a zero pivot, the dense kernel also scales them by its infinite reciprocal (faer-010). + +## How the fork fixed it + +The fork adds `LuError::NumericSingular { index }`. This is the variant the upstream maintainer asked for in upstream PR #332, which was closed unmerged. +- The simplicial backend returns it when the chosen pivot is exactly zero. +- The supernodal backend returns it when a diagonal entry of a factored panel is exactly zero. + +`index` counts elimination steps in the fill-reducing column order, as `SymbolicSingular`'s does. Tolerance-based detection is not part of the fix; it is tracked separately. Run against the fork, both backends report `NumericSingular`. + +Fork regression tests: `sparse::linalg::lu::tests::codeberg_339` (`f32`, `f64`, `c32`, `c64`, both backends) and `sparse::linalg::lu::tests::numeric_singular_after_an_elimination_step`, and `sparse::linalg::lu::tests::nonzero_subnormal_pivot_is_not_singular` (a nonzero subnormal complex pivot is not reported). + +## 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 + +* Upstream report: sarah-quinones/faer#339 (simplicial only). Upstream PR: sarah-quinones/faer#332. +* Fork fix: tensor4all/faer#132. diff --git a/upstream-reports/faer-028-sparse-lu-numeric-singular/src/main.rs b/upstream-reports/faer-028-sparse-lu-numeric-singular/src/main.rs new file mode 100644 index 00000000..b07a9b09 --- /dev/null +++ b/upstream-reports/faer-028-sparse-lu-numeric-singular/src/main.rs @@ -0,0 +1,66 @@ +// Sparse LU of the singular matrix A = [[1, 1], [1, 1]], with each backend forced. +// The numeric factorization must report the zero pivot instead of returning Ok. +use dyn_stack::{MemBuffer, MemStack}; +use faer::prelude::*; +use faer::sparse::linalg::lu::{factorize_symbolic_lu, LuSymbolicParams, NumericLu}; +use faer::sparse::linalg::SupernodalThreshold; +use faer::sparse::{SparseColMat, Triplet}; +use faer::Conj; + +fn main() { + let entries = [ + Triplet::new(0, 0, 1.0f64), + Triplet::new(0, 1, 1.0), + Triplet::new(1, 0, 1.0), + Triplet::new(1, 1, 1.0), + ]; + let a = SparseColMat::::try_new_from_triplets(2, 2, &entries).unwrap(); + let mut wrong = false; + for (name, threshold) in [ + ("simplicial", SupernodalThreshold::FORCE_SIMPLICIAL), + ("supernodal", SupernodalThreshold::FORCE_SUPERNODAL), + ] { + let symbolic = factorize_symbolic_lu( + a.symbolic(), + LuSymbolicParams { + supernodal_flop_ratio_threshold: threshold, + ..Default::default() + }, + ) + .unwrap(); + let mut numeric = NumericLu::::new(); + let result = symbolic.factorize_numeric_lu( + &mut numeric, + a.as_ref(), + Par::Seq, + MemStack::new(&mut MemBuffer::new( + symbolic.factorize_numeric_lu_scratch::(Par::Seq, Default::default()), + )), + Default::default(), + ); + match result { + Err(e) => println!("{name}: expected an error, got {e:?}"), + Ok(lu) => { + let mut x = Mat::::ones(2, 1); + lu.solve_in_place_with_conj( + Conj::No, + x.as_mut(), + Par::Seq, + MemStack::new(&mut MemBuffer::new( + symbolic.solve_in_place_scratch::(1, Par::Seq), + )), + ); + println!( + "{name}: expected an error, got Ok; solve of A x = [1, 1] gives x = [{}, {}]", + x[(0, 0)], + x[(1, 0)] + ); + wrong = true; + } + } + } + if wrong { + println!("WRONG: a singular matrix factors without an error"); + std::process::exit(1); + } +}