operator: final residual on NoConvergence, BiCGSTAB breakdown stops at the last finite iterate (#133, #134)

PATCHSET.md

@@ -106,6 +106,7 @@ gate has changed and two kinds of adoption are now in the ledger.
 | 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). |
 | 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. |
 | 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/operator/bicgstab.rs

@@ -149,6 +149,11 @@ pub fn bicgstab<T: ComplexField>(
 				non_exhaustive: NonExhaustive(()),
 			});
 		}
+		if params.initial_guess == InitialGuessStatus::Zero {
+			// every exit returns `x`, so an implicit zero is written out once
+			// here
+			x.fill(zero());
+		}
 		let rel_threshold = params.rel_tolerance * b_norm;
 		let abs_threshold = params.abs_tolerance;
 		let threshold = if abs_threshold > rel_threshold {
@@ -170,7 +175,7 @@ pub fn bicgstab<T: ComplexField>(
 		let (mut r_tilde, stack) =
 			unsafe { temp_mat_uninit::<T, _, _>(n, k, stack) };
 		let mut r_tilde = r_tilde.as_mat_mut();
-		let abs_residual =
+		let mut abs_residual =
 			if params.initial_guess == InitialGuessStatus::MaybeNonZero {
 				A.apply(r.rb_mut(), x.rb(), par.reborrow(), stack);
 				z!(&mut r, &b).for_each(|uz!(r, b)| *r = b - &*r);
@@ -223,16 +228,36 @@ pub fn bicgstab<T: ComplexField>(
 					stack,
 					Default::default(),
 				);
+			// the tolerance is relative to what an entry of `r̃ᵀ v` could be,
+			// not to its own largest pivot: otherwise a zero first pivot
+			// passes and the step divides by zero
+			let max_col_norm = |m: MatRef<'_, T>| {
+				let mut max = zero::<T::Real>();
+				for j in 0..k {
+					let norm = m.col(j).norm_l2();
+					if norm > max {
+						max = norm;
+					}
+				}
+				max
+			};
 			let mut rank = k;
 			let tol = eps::<T::Real>()
 				* from_f64::<T::Real>(k as f64)
-				* rtv[(0, 0)].abs();
+				* max_col_norm(r_tilde.rb())
+				* max_col_norm(v.rb());
 			for i in 0..k {
-				if rtv[(i, i)].abs() < tol {
+				if !(rtv[(i, i)].abs() > tol) {
 					rank = i;
 					break;
 				}
 			}
+			if rank == 0 {
+				// breakdown: no direction left; `x` is the last iterate and `r`
+				// its residual
+				break;
+			}
+			let s_norm;
 			let (mut s, stack) =
 				unsafe { temp_mat_uninit::<T, _, _>(n, k, stack) };
 			let mut s = s.as_mat_mut();
@@ -278,6 +303,11 @@ pub fn bicgstab<T: ComplexField>(
 					-one::<T>(),
 					par.reborrow(),
 				);
+				s_norm = s.norm_l2();
+				if !s_norm.is_finite() {
+					// the step is not usable: keep `x` and `r` as they are
+					break;
+				}
 				crate::linalg::matmul::matmul(
 					x.rb_mut(),
 					if iter == 0
@@ -293,11 +323,10 @@ pub fn bicgstab<T: ComplexField>(
 					par.reborrow(),
 				);
 			}
-			let norm = s.norm_l2();
-			if norm < threshold {
+			if s_norm < threshold {
 				return Ok(BicgInfo {
-					rel_residual: &norm / b_norm,
-					abs_residual: norm,
+					rel_residual: &s_norm / b_norm,
+					abs_residual: s_norm,
 					iter_count: iter + 1,
 					non_exhaustive: NonExhaustive(()),
 				});
@@ -337,10 +366,16 @@ pub fn bicgstab<T: ComplexField>(
 				compute_w(kt.rb(), ks.rb())
 			};
 			z!(&mut r, &s, &t).for_each(|uz!(r, s, t)| *r = s - w * t);
+			let norm = r.norm_l2();
+			if !norm.is_finite() {
+				// `x` keeps the first half step, whose residual is `s`
+				abs_residual = s_norm;
+				break;
+			}
 			z!(&mut x, &z).for_each(|uz!(x, z)| *x += w * z);
 			z!(&mut p, &v).for_each(|uz!(p, v)| *p -= w * v);
 			callback(x.rb());
-			let norm = r.norm_l2();
+			abs_residual = norm.copy();
 			if norm < threshold {
 				return Ok(BicgInfo {
 					rel_residual: &norm / b_norm,
@@ -487,4 +522,246 @@ mod tests {
 			assert!((out[(0, 0)] - 1.0).abs() < 1e-10);
 		}
 	}
+	#[test]
+	fn codeberg_328() {
+		// one step does not converge: the error carries the residual of the
+		// returned `x`
+		let ref A = mat![[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 3.0]];
+		let ref rhs = mat![[1.0], [1.0], [1.0]];
+		let precond = IdentityPrecond { dim: 3 };
+		let params = BicgParams {
+			max_iters: 1,
+			..BicgParams::default()
+		};
+		let ref mut out = Mat::<f64>::zeros(3, 1);
+		let result = bicgstab(
+			out.as_mut(),
+			precond,
+			precond,
+			A.as_ref(),
+			rhs.as_ref(),
+			params,
+			|_| {},
+			Par::Seq,
+			MemStack::new(&mut MemBuffer::new(bicgstab_scratch(
+				precond,
+				precond,
+				A.as_ref(),
+				1,
+				ParShape::Seq,
+			))),
+		);
+		let BicgError::NoConvergence {
+			abs_residual,
+			rel_residual,
+		} = result.unwrap_err();
+		let true_residual = (rhs - A * &*out).norm_l2();
+		assert!((abs_residual - true_residual).abs() < 1e-12);
+		assert!((rel_residual - true_residual / rhs.norm_l2()).abs() < 1e-12);
+		assert!(abs_residual < 0.5);
+	}
+	#[test]
+	fn codeberg_327() {
+		// after one step `r̃ᵀ r = 0`, so the next `r̃ᵀ v` is zero: a breakdown
+		let ref A = mat![[1.0, 0.0, 0.0], [-1.0, 1.0, 0.0], [0.0, -1.0, 1.0]];
+		let ref rhs = mat![[1.0], [0.0], [0.0]];
+		let precond = IdentityPrecond { dim: 3 };
+		for initial_guess in
+			[InitialGuessStatus::MaybeNonZero, InitialGuessStatus::Zero]
+		{
+			let params = BicgParams {
+				initial_guess,
+				max_iters: 10,
+				..BicgParams::default()
+			};
+			let ref mut out = Mat::<f64>::zeros(3, 1);
+			let result = bicgstab(
+				out.as_mut(),
+				precond,
+				precond,
+				A.as_ref(),
+				rhs.as_ref(),
+				params,
+				|_| {},
+				Par::Seq,
+				MemStack::new(&mut MemBuffer::new(bicgstab_scratch(
+					precond,
+					precond,
+					A.as_ref(),
+					1,
+					ParShape::Seq,
+				))),
+			);
+			let BicgError::NoConvergence { abs_residual, .. } =
+				result.unwrap_err();
+			assert!(out.norm_l2().is_finite());
+			let true_residual = (rhs - A * &*out).norm_l2();
+			assert!((abs_residual - true_residual).abs() < 1e-12);
+		}
+	}
+	#[test]
+	fn breakdown_when_the_operator_maps_the_direction_to_zero() {
+		// `A r̃ = 0`, so `v = 0` and the tolerance is zero: still a breakdown
+		let ref A = mat![[0.0, 0.0], [0.0, 1.0]];
+		let ref rhs = mat![[1.0], [0.0]];
+		let precond = IdentityPrecond { dim: 2 };
+		let params = BicgParams {
+			max_iters: 10,
+			..BicgParams::default()
+		};
+		let ref mut out = Mat::<f64>::zeros(2, 1);
+		let result = bicgstab(
+			out.as_mut(),
+			precond,
+			precond,
+			A.as_ref(),
+			rhs.as_ref(),
+			params,
+			|_| {},
+			Par::Seq,
+			MemStack::new(&mut MemBuffer::new(bicgstab_scratch(
+				precond,
+				precond,
+				A.as_ref(),
+				1,
+				ParShape::Seq,
+			))),
+		);
+		let BicgError::NoConvergence { abs_residual, .. } = result.unwrap_err();
+		assert!(out.norm_l2() == 0.0);
+		assert!(abs_residual == 1.0);
+	}
+	#[test]
+	fn breakdown_on_the_first_step_with_a_zero_guess() {
+		// `r̃ᵀ A r̃ = 0` at once; a zero initial guess must come back as zeros,
+		// not as the caller's buffer
+		let ref A = mat![[0.0, 1.0], [-1.0, 0.0]];
+		let ref rhs = mat![[1.0], [0.0]];
+		let precond = IdentityPrecond { dim: 2 };
+		let params = BicgParams {
+			initial_guess: InitialGuessStatus::Zero,
+			max_iters: 10,
+			..BicgParams::default()
+		};
+		let ref mut out = Mat::<f64>::full(2, 1, f64::NAN);
+		let result = bicgstab(
+			out.as_mut(),
+			precond,
+			precond,
+			A.as_ref(),
+			rhs.as_ref(),
+			params,
+			|_| {},
+			Par::Seq,
+			MemStack::new(&mut MemBuffer::new(bicgstab_scratch(
+				precond,
+				precond,
+				A.as_ref(),
+				1,
+				ParShape::Seq,
+			))),
+		);
+		let BicgError::NoConvergence { abs_residual, .. } = result.unwrap_err();
+		assert!(out.norm_l2() == 0.0);
+		assert!(abs_residual == 1.0);
+	}
+	#[test]
+	fn nonfinite_second_half_step_keeps_the_first() {
+		// `A s = 0`, so `w` is 0 / 0: `x` keeps the alpha step, whose residual
+		// is `s`
+		let ref A = mat![[2.0, 0.0], [1.0, 0.0]];
+		let ref rhs = mat![[1.0], [0.0]];
+		let precond = IdentityPrecond { dim: 2 };
+		let params = BicgParams {
+			max_iters: 10,
+			..BicgParams::default()
+		};
+		let ref mut out = Mat::<f64>::zeros(2, 1);
+		let result = bicgstab(
+			out.as_mut(),
+			precond,
+			precond,
+			A.as_ref(),
+			rhs.as_ref(),
+			params,
+			|_| {},
+			Par::Seq,
+			MemStack::new(&mut MemBuffer::new(bicgstab_scratch(
+				precond,
+				precond,
+				A.as_ref(),
+				1,
+				ParShape::Seq,
+			))),
+		);
+		let BicgError::NoConvergence { abs_residual, .. } = result.unwrap_err();
+		assert!(out.norm_l2().is_finite());
+		assert!((abs_residual - (rhs - A * &*out).norm_l2()).abs() < 1e-12);
+		assert!((abs_residual - 0.5).abs() < 1e-12);
+	}
+	#[test]
+	fn zero_guess_is_written_out_without_an_iteration() {
+		let ref A = mat![[1.0]];
+		let ref rhs = mat![[1.0]];
+		let precond = IdentityPrecond { dim: 1 };
+		let params = BicgParams {
+			initial_guess: InitialGuessStatus::Zero,
+			max_iters: 0,
+			..BicgParams::default()
+		};
+		let ref mut out = Mat::<f64>::full(1, 1, 7.0);
+		let result = bicgstab(
+			out.as_mut(),
+			precond,
+			precond,
+			A.as_ref(),
+			rhs.as_ref(),
+			params,
+			|_| {},
+			Par::Seq,
+			MemStack::new(&mut MemBuffer::new(bicgstab_scratch(
+				precond,
+				precond,
+				A.as_ref(),
+				1,
+				ParShape::Seq,
+			))),
+		);
+		let BicgError::NoConvergence { abs_residual, .. } = result.unwrap_err();
+		assert!(out[(0, 0)] == 0.0);
+		assert!(abs_residual == 1.0);
+	}
+	#[test]
+	fn block_rhs_with_identity_is_not_a_breakdown() {
+		// `r̃ᵀ v` is the identity: its pivots are 1, whatever the block width
+		let k = 64;
+		let ref A = Mat::<f32>::identity(k, k);
+		let ref rhs = Mat::<f32>::identity(k, k);
+		let precond = IdentityPrecond { dim: k };
+		let params = BicgParams {
+			initial_guess: InitialGuessStatus::Zero,
+			max_iters: 1,
+			..BicgParams::default()
+		};
+		let ref mut out = Mat::<f32>::zeros(k, k);
+		let result = bicgstab(
+			out.as_mut(),
+			precond,
+			precond,
+			A.as_ref(),
+			rhs.as_ref(),
+			params,
+			|_| {},
+			Par::Seq,
+			MemStack::new(&mut MemBuffer::new(bicgstab_scratch(
+				precond,
+				precond,
+				A.as_ref(),
+				k,
+				ParShape::Seq,
+			))),
+		);
+		assert!(result.is_ok());
+		assert!((A * &*out - rhs).norm_l2() < 1e-5);
+	}
 }