sparse lu: report an exact zero pivot as LuError::NumericSingular in both backends (#132)

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/). |

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::<T::Real>() {
+					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::<T::Real>() {
+				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<T: ComplexField>(
+		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::<usize, T>::new();
+		symbolic
+			.factorize_numeric_lu(
+				&mut numeric,
+				A,
+				Par::Seq,
+				MemStack::new(&mut MemBuffer::new(
+					symbolic.factorize_numeric_lu_scratch::<T>(
+						ParShape::Seq,
+						Default::default(),
+					),
+				)),
+				Default::default(),
+			)
+			.map(|_| ())
+	}
+	fn numeric_singular_index<T: ComplexField>(
+		n: usize,
+		entries: &[(usize, usize, f64)],
+		supernodal_flop_ratio_threshold: SupernodalThreshold,
+	) -> Option<usize> {
+		let triplets = entries
+			.iter()
+			.map(|&(i, j, v)| Triplet::new(i, j, from_f64::<T>(v)))
+			.collect::<alloc::vec::Vec<_>>();
+		let A =
+			SparseColMat::<usize, T>::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::<f64>(2, &A, threshold) == Some(1));
+			assert!(numeric_singular_index::<f32>(2, &A, threshold) == Some(1));
+			assert!(numeric_singular_index::<c64>(2, &A, threshold) == Some(1));
+			assert!(numeric_singular_index::<c32>(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::<f64>(4, &A, threshold),
+				Some(index) if index < 4
+			));
+			assert!(numeric_singular_index::<f64>(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::<usize, c64>::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());
+		}
+	}
 }

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),
 }

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).