evd: the 2 x 2 block solves of the real eigenvector back substitution follow ?LALN2

faer/src/linalg/evd/mod.rs

@@ -901,14 +901,19 @@ fn solve_real_shifted_upper_quasi_triangular_system<T: RealField>(
 				);
 				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<T: RealField>(
 		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<T: RealField>(
+	a: (T, T),
+	b: T,
+	c: T,
+	r0: (T, T),
+	r1: (T, T),
+	smin: &T,
+) -> ((T, T), (T, T)) {
+	type C<T> = (T, T);
+	fn abs1<T: RealField>(z: &C<T>) -> T {
+		z.0.abs() + z.1.abs()
+	}
+	fn mul<T: RealField>(x: &C<T>, y: &C<T>) -> C<T> {
+		(&x.0 * &y.0 - &x.1 * &y.1, &x.0 * &y.1 + &x.1 * &y.0)
+	}
+	fn sub<T: RealField>(x: &C<T>, y: &C<T>) -> C<T> {
+		(&x.0 - &y.0, &x.1 - &y.1)
+	}
+	fn div<T: RealField>(x: &C<T>, y: &C<T>) -> C<T> {
+		// 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::<T>;
+	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<T: RealField>(
 	A: MatRef<'_, T>,
 	p: T,
@@ -1043,25 +1112,24 @@ fn solve_cplx_shifted_upper_quasi_triangular_system<T: RealField>(
 				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::<T>(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 {

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::<Scalar>::zeros(n, n);
+		let mut r = Mat::<Scalar>::zeros(n, n);
+		let mut wr = Col::<Real>::zeros(n);
+		let mut wi = Col::<Real>::zeros(n);
+		let p = Spec::new(EvdParams { ..auto!(Scalar) });
+		let mut mem = MemBuffer::new(evd_scratch::<Scalar>(
+			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 !<Scalar as ComplexField>::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::<Scalar>::zeros(4, 4);
+	let p = Spec::new(EvdFromSchurParams { ..auto!(Scalar) });
+	let mut mem = MemBuffer::new(eigenvectors_from_schur_scratch::<Scalar>(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 !<Scalar as ComplexField>::IS_REAL {
+		return;
+	}
+	let t = if core::mem::size_of::<Real>() == 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::<Scalar>::zeros(3, 3);
+	let p = Spec::new(EvdFromSchurParams { ..auto!(Scalar) });
+	let mut mem = MemBuffer::new(eigenvectors_from_schur_scratch::<Scalar>(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}");
+}

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