From f27f714ab9e8d6d47604f576d105a86b7a846713 Mon Sep 17 00:00:00 2001
From: Kristof Schroeder <kristof_schroeder@web.de>
Date: Wed, 23 Sep 2026 14:27:31 +0200
Subject: [PATCH 1/4] =?UTF-8?q?fix(lsmr):=20run=20the=20recurrence=20on=20?=
 =?UTF-8?q?u=E2=82=81=20and=20test=20its=20stop=20in=20exact=20exponents?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

At ‖A‖‖b‖ > 1e308 the unscaled `ζ̄₁ = α₁β₁` overflowed, and `‖A‖²` in
the Frobenius estimate left the range from ‖A‖ ≈ 1e±154. The scalar
recurrence now runs on `u₁ = b / β₁`, with β₁ restored only in the
solution update and the reported residual; ‖A‖_F accumulates by hypot.

The backward error `‖Aᵀr‖ / (‖A‖‖r‖)` and the reported relative
normal-equation residual are formed as `Magnitude`s (m·2^e), whose
products and ratios keep their exponent exactly, so neither overflows
where the ratio itself is representable.

A power-of-two scaling sweep pins the unpreconditioned solve to its
unit-scale twin at ‖A‖ = 2^±5xx and ‖A‖‖b‖ = 2^±1030.
---
 CHANGELOG.md                                  |   2 +-
 crates/schwarz-precond/src/lsmr.rs            |  36 ++--
 crates/schwarz-precond/src/lsmr/magnitude.rs  | 159 ++++++++++++++++++
 crates/schwarz-precond/src/lsmr/recurrence.rs |  75 ++++-----
 .../schwarz-precond/src/lsmr/tests/audit.rs   |   9 +-
 .../schwarz-precond/src/lsmr/tests/solve.rs   |  66 ++++++--
 6 files changed, 281 insertions(+), 66 deletions(-)
 create mode 100644 crates/schwarz-precond/src/lsmr/magnitude.rs

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 46c66497..a7270eb5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -43,7 +43,7 @@ and this project follows [Semantic Versioning](https://semver.org/).
 - LSMR no longer certifies a stop whose recurrence estimates collapsed. Tolerance stops are checked against `‖b − A x‖`, a non-finite `α`, `β`, `⟨v, Mv⟩`, or `‖b‖` fails with `SolveError::InvalidInput`, and an overflowing or subnormal `‖A‖` no longer zeroes the normal-equation ratio; a failed check reports `LsmrStopReason::FalseConvergence` with `converged = false` and the returned iterate's normal-equation residual (#290, #297, #303, #362).
 - A warm-started solve measures its residuals against the original `b`, including at a budget stop and where `‖Aᵀb‖²` overflows.
 - An `α` or `β` whose square underflows is recovered from a scaled norm, where LSMR reported `x = 0` converged.
-- LSMR's solution update no longer drops out when `‖A‖` is small (`≈ 1e-9`) or beyond `1e±154`, where it reported `x = 0` converged.
+- LSMR's solution update no longer drops out when `‖A‖` is small (`≈ 1e-9`) or beyond `1e±154`, where it reported `x = 0` converged, and a solve whose `‖A‖²` or `‖A‖‖b‖` leaves the double range no longer fails with `SolveError::InvalidInput`, misses its normal-equation stop, or reports a zero normal-equation residual.
 - A preconditioned solve with local reorthogonalization (`local_size`) could fail with `SolveError::InvalidInput` ("preconditioner not positive definite") on a positive definite preconditioner, once reorthogonalization cancelled the last Krylov direction to rounding noise; the pair is now recomputed before the check.
 - LSMR's residual estimate is its own `‖r_k‖` rather than LSQR's smaller `|φ̄_k|`, which let `ResidualTolerance` fire before the tolerance was met.
 - A non-finite or negative `ScalingConfig::tolerance` silently disabled the dominance certificate under both failure policies, since every comparison against it is `>`; it is now rejected as `BuildError::InvalidScalingTolerance`.
diff --git a/crates/schwarz-precond/src/lsmr.rs b/crates/schwarz-precond/src/lsmr.rs
index 87dcb0fe..532e961a 100644
--- a/crates/schwarz-precond/src/lsmr.rs
+++ b/crates/schwarz-precond/src/lsmr.rs
@@ -7,6 +7,7 @@
 mod bidiag;
 #[cfg(test)]
 mod fixtures;
+mod magnitude;
 mod recurrence;
 #[cfg(test)]
 mod tests;
@@ -18,6 +19,7 @@ use bidiag::{
     axpby, metric_gradient_norm, residual_into, BidiagStep, Bidiagonalization, GolubKahan,
     ModifiedGolubKahan,
 };
+use magnitude::Magnitude;
 use recurrence::{ConvergenceCriteria, LsmrRecurrenceState, RotationStep, SolutionState, Stop};
 
 /// Euclidean norm of a vector.
@@ -333,7 +335,7 @@ pub fn mlsmr<A: Operator + ?Sized, M: Operator + ?Sized>(
         ModifiedGolubKahan::init(operator, preconditioner, &rhs, rhs_norm, local_size)?;
     let warm_start = warm_start.zip(metric).map(|(x0, metric)| WarmStart {
         x0,
-        reference: NormalEqReference::warm(metric, step1),
+        reference: NormalEqReference::warm(Magnitude::from(metric), step1),
     });
     let reference_norm = if b_norm > 0.0 { b_norm } else { rhs_norm };
     let criteria = ConvergenceCriteria::new(reference_norm, tol);
@@ -345,25 +347,31 @@ const MAX_RESTARTS: usize = 2;
 
 /// `‖Aᵀb‖` in the stream's metric, fixed for the solve so every pass reports against it.
 #[derive(Clone, Copy)]
-struct NormalEqReference(f64);
+struct NormalEqReference(Magnitude);
 
 impl NormalEqReference {
-    /// A cold stream's first `(α₁, β₁)` is `‖Aᵀb‖` already.
+    /// A cold stream's first `(α₁, β₁)` is `‖Aᵀb‖` already, though the product may not be a double.
     fn cold(step1: BidiagStep) -> Self {
-        Self((step1.alpha * step1.beta).abs().max(f64::MIN_POSITIVE))
+        Self(Magnitude::product(step1.alpha, step1.beta))
     }
 
     /// A reference that carries no information leaves the stream's own `ζ̄₀` to divide.
-    fn warm(metric: f64, step1: BidiagStep) -> Self {
-        if metric > 0.0 && metric.is_finite() {
+    fn warm(metric: Magnitude, step1: BidiagStep) -> Self {
+        if metric.is_normal() {
             Self(metric)
         } else {
             Self::cold(step1)
         }
     }
 
-    fn relative(self, estimate: f64) -> f64 {
-        estimate / self.0
+    /// A vanished reference divides as the smallest normal instead.
+    fn relative(self, estimate: Magnitude) -> f64 {
+        let reference = if self.0 > Magnitude::ZERO {
+            self.0
+        } else {
+            Magnitude::from(f64::MIN_POSITIVE)
+        };
+        (estimate / reference).to_f64()
     }
 }
 
@@ -397,13 +405,17 @@ fn lsmr_from_bidiag<B: Bidiagonalization>(
             let (converged, normal_eq_residual) = match bidiag.hidden_gradient()? {
                 None => (true, 0.0),
                 Some(per_unit_residual) => {
-                    let normar = step1.beta * per_unit_residual;
+                    let normar = Magnitude::product(step1.beta, per_unit_residual);
                     let plain = bidiag.plain_gradient(b)?;
                     let a_norm_below = plain / vec_norm(b).max(f64::MIN_POSITIVE);
                     let informative = plain > 0.0 && plain.is_finite();
                     (
                         criteria.corroborates(step1.beta, normar, a_norm_below),
-                        if informative { normar / plain } else { 1.0 },
+                        if informative {
+                            (normar / Magnitude::from(plain)).to_f64()
+                        } else {
+                            1.0
+                        },
                     )
                 }
             };
@@ -426,7 +438,7 @@ fn lsmr_from_bidiag<B: Bidiagonalization>(
         let mut escalation = escalation.map(EscalationPolicy::handler);
         let mut convergence = criteria.start(step1.alpha);
         let mut recurrence = LsmrRecurrenceState::init(step1);
-        let mut solution = SolutionState::init(bidiag.v());
+        let mut solution = SolutionState::init(bidiag.v(), step1.beta);
         let mut prev_rot = RotationStep::initial();
         let stop_reason = 'pass: {
             while iterations < maxiter {
@@ -493,7 +505,7 @@ fn lsmr_from_bidiag<B: Bidiagonalization>(
                 // The estimate is the refuted claim; a seed from the staged residual measures `x`.
                 result.normal_eq_residual = if residual_norm.is_finite() {
                     let step = bidiag.restart(residual_norm)?;
-                    reference.relative(step.alpha * step.beta)
+                    reference.relative(Magnitude::product(step.alpha, step.beta))
                 } else {
                     residual_norm
                 };
diff --git a/crates/schwarz-precond/src/lsmr/magnitude.rs b/crates/schwarz-precond/src/lsmr/magnitude.rs
new file mode 100644
index 00000000..0e393480
--- /dev/null
+++ b/crates/schwarz-precond/src/lsmr/magnitude.rs
@@ -0,0 +1,159 @@
+//! Norms as `m · 2^e`, for comparing gradients of problems whose scales share no range.
+
+use std::cmp::Ordering;
+use std::ops::{Div, Mul};
+
+/// A nonnegative `m · 2^e`, `m ∈ [1, 2)`, zero or non-finite, for expressions a few factors deep.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub(super) struct Magnitude {
+    m: f64,
+    e: i32,
+}
+
+impl Magnitude {
+    pub(super) const ZERO: Self = Self { m: 0.0, e: 0 };
+
+    /// `a · b`, exact in exponent where the `f64` product would over- or underflow.
+    pub(super) fn product(a: f64, b: f64) -> Self {
+        Self::from(a) * Self::from(b)
+    }
+
+    /// The nearest `f64`, saturating to 0 or ∞; a subnormal result rounds a second time.
+    pub(super) fn to_f64(self) -> f64 {
+        if !self.m.is_normal() || self.e >= f64::MAX_EXP {
+            return if self.m.is_normal() {
+                f64::INFINITY
+            } else {
+                self.m
+            };
+        }
+        if self.e >= f64::MIN_EXP - 1 {
+            return self.m * pow2(self.e);
+        }
+        // A subnormal result: the first factor is exact, so only the second rounds.
+        self.m * pow2((self.e + 64).max(f64::MIN_EXP - 1)) * pow2(-64)
+    }
+
+    /// Neither zero, infinite nor NaN.
+    pub(super) fn is_normal(self) -> bool {
+        self.m.is_normal()
+    }
+
+    fn normalized(m: f64, e: i32) -> Self {
+        if !m.is_normal() {
+            return Self { m, e: 0 };
+        }
+        let Self { m, e: shift } = Self::from(m);
+        Self { m, e: e + shift }
+    }
+}
+
+impl From<f64> for Magnitude {
+    fn from(x: f64) -> Self {
+        if !x.is_finite() || x == 0.0 {
+            return Self { m: x, e: 0 };
+        }
+        // A subnormal carries its leading bit below the mantissa field; 2^64 lifts it exactly.
+        let (x, bias) = if x.is_normal() {
+            (x, 0)
+        } else {
+            (x * pow2(64), -64)
+        };
+        debug_assert!(x > 0.0, "a norm, got {x:e}");
+        let bits = x.to_bits();
+        let e = ((bits >> 52) & 0x7ff) as i32 - 1023;
+        let m = f64::from_bits((bits & ((1 << 52) - 1)) | (1023 << 52));
+        Self { m, e: e + bias }
+    }
+}
+
+impl Mul for Magnitude {
+    type Output = Self;
+    fn mul(self, rhs: Self) -> Self {
+        Self::normalized(self.m * rhs.m, self.e + rhs.e)
+    }
+}
+
+impl Div for Magnitude {
+    type Output = Self;
+    fn div(self, rhs: Self) -> Self {
+        Self::normalized(self.m / rhs.m, self.e - rhs.e)
+    }
+}
+
+impl PartialOrd for Magnitude {
+    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
+        if self.m.is_normal() && other.m.is_normal() {
+            return Some(self.e.cmp(&other.e).then(self.m.total_cmp(&other.m)));
+        }
+        // 0, ∞ and NaN order against any `m ∈ [1, 2)` as they do against its value.
+        self.m.partial_cmp(&other.m)
+    }
+}
+
+/// `2^e` for a normal exponent, built from its bits since `powi` underflows on the way there.
+fn pow2(e: i32) -> f64 {
+    debug_assert!((f64::MIN_EXP - 1..f64::MAX_EXP).contains(&e));
+    f64::from_bits(((e + 1023) as u64) << 52)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::Magnitude;
+    use rstest::rstest;
+
+    #[rstest]
+    #[case::one(1.0)]
+    #[case::large(3.7e300)]
+    #[case::smallest_normal(f64::MIN_POSITIVE)]
+    #[case::subnormal(1.3e-310)]
+    #[case::smallest_subnormal(f64::from_bits(1))]
+    #[case::largest(f64::MAX)]
+    #[case::zero(0.0)]
+    #[case::infinity(f64::INFINITY)]
+    fn a_float_round_trips(#[case] x: f64) {
+        assert_eq!(Magnitude::from(x).to_f64().to_bits(), x.to_bits());
+    }
+
+    /// Where the `f64` arithmetic stays in range the two agree bitwise.
+    #[rstest]
+    #[case::unit(3.0, 7.0, 5.0)]
+    #[case::large(3e150, 7e150, 5e-7)]
+    #[case::small(3e-150, 7e-150, 5e7)]
+    fn in_range_arithmetic_is_the_float_arithmetic(#[case] a: f64, #[case] b: f64, #[case] c: f64) {
+        let exact = Magnitude::product(a, b) / Magnitude::from(c);
+        assert_eq!(exact.to_f64().to_bits(), (a * b / c).to_bits());
+    }
+
+    /// A quotient back in range survives a product that is not.
+    #[rstest]
+    #[case::overflowing_product(1e300, 3e300, 1e300, 3e300)]
+    #[case::underflowing_product(1e-300, 3e-300, 1e-300, 3e-300)]
+    #[case::subnormal_quotient(1e-300, 3e-10, 1.0, 3e-310)]
+    fn out_of_range_products_divide_back(
+        #[case] a: f64,
+        #[case] b: f64,
+        #[case] c: f64,
+        #[case] expected: f64,
+    ) {
+        let got = (Magnitude::product(a, b) / Magnitude::from(c)).to_f64();
+        assert!(
+            (got / expected - 1.0).abs() < 1e-9,
+            "{got:e} vs {expected:e}"
+        );
+    }
+
+    #[test]
+    fn the_order_is_the_value_order() {
+        let tiny = Magnitude::product(1e-300, 1e-300);
+        let huge = Magnitude::product(1e300, 1e300);
+        assert!(Magnitude::ZERO < tiny && tiny < Magnitude::from(1.0));
+        assert!(Magnitude::from(1.0) < huge && huge < Magnitude::from(f64::INFINITY));
+        assert!(Magnitude::from(1.5) > Magnitude::from(1.25));
+        assert_eq!(tiny.to_f64(), 0.0);
+        assert_eq!(huge.to_f64(), f64::INFINITY);
+        assert!(Magnitude::from(f64::NAN)
+            .partial_cmp(&Magnitude::ZERO)
+            .is_none());
+    }
+}
diff --git a/crates/schwarz-precond/src/lsmr/recurrence.rs b/crates/schwarz-precond/src/lsmr/recurrence.rs
index e58c8619..150ffec7 100644
--- a/crates/schwarz-precond/src/lsmr/recurrence.rs
+++ b/crates/schwarz-precond/src/lsmr/recurrence.rs
@@ -1,11 +1,12 @@
 //! LSMR scalar/vector recurrence consuming the bidiagonalization stream.
 //!
 //! Given `(α, β)` pairs from a [`super::bidiag::Bidiagonalization`], this
-//! module builds the two interleaved Givens rotation chains (P̂_k, P̄_k)
+//! module builds the three Givens rotation chains (P̂_k, P̄_k, Q̃_k)
 //! that yield Algorithm 2.8 of Fong & Saunders, advances the `(x, h, h̄)`
 //! solution recurrence, and tracks the dual stopping criterion.
 
 use super::bidiag::{BidiagStep, LSMR_PAR_THRESHOLD, LSMR_UPDATE_CHUNK};
+use super::magnitude::Magnitude;
 use rayon::iter::{IndexedParallelIterator, ParallelIterator};
 use rayon::prelude::{ParallelSlice, ParallelSliceMut};
 
@@ -54,8 +55,10 @@ impl RotationStep {
     }
 }
 
-/// LSMR scalar state: the three rotation sequences of Fong & Saunders.
+/// LSMR scalar state: the three rotation sequences of Fong & Saunders, run on `u₁ = b / β₁`.
 pub(super) struct LsmrRecurrenceState {
+    /// `β₁`; on `b` itself `ζ̄₁ = α₁β₁` leaves the double range while `α₁` and `β₁` do not.
+    rhs_norm: f64,
     bidiag_qr: BidiagQr,
     normal_eq_qr: NormalEqQr,
     residual: ResidualChain,
@@ -64,15 +67,16 @@ pub(super) struct LsmrRecurrenceState {
 impl LsmrRecurrenceState {
     pub(super) fn init(s1: BidiagStep) -> Self {
         Self {
+            rhs_norm: s1.beta,
             bidiag_qr: BidiagQr {
                 alpha_bar: s1.alpha,
-                beta_dd: s1.beta,
+                beta_dd: 1.0,
             },
             normal_eq_qr: NormalEqQr {
                 c_bar: 1.0,
                 s_bar: 0.0,
-                zeta_bar: s1.alpha * s1.beta,
-                zeta0: (s1.alpha * s1.beta).abs().max(f64::MIN_POSITIVE),
+                zeta_bar: s1.alpha,
+                zeta0: s1.alpha,
             },
             residual: ResidualChain {
                 beta_d: 0.0,
@@ -94,15 +98,21 @@ impl LsmrRecurrenceState {
 
     /// Running estimate of LSMR's own `‖r_k‖`.
     pub(super) fn residual_estimate(&self) -> f64 {
+        self.rhs_norm * self.unit_residual()
+    }
+
+    /// `‖r_k‖` for `u₁`, where `‖r_k‖ ≤ 1`.
+    fn unit_residual(&self) -> f64 {
         self.residual.normr(self.bidiag_qr.beta_dd)
     }
 
-    pub(super) fn normal_eq_residual_estimate(&self) -> f64 {
-        self.normal_eq_qr.normar()
+    /// Running estimate of `‖Âᵀ r_k‖` for `b`, a product that may leave the double range.
+    pub(super) fn normal_eq_residual_estimate(&self) -> Magnitude {
+        Magnitude::product(self.rhs_norm, self.normal_eq_qr.normar())
     }
 
     pub(super) fn relative_normal_eq_residual(&self) -> f64 {
-        self.normal_eq_qr.relative_normar()
+        self.normal_eq_qr.normar() / self.normal_eq_qr.zeta0
     }
 }
 
@@ -141,7 +151,7 @@ struct NormalEqQr {
     c_bar: f64,
     s_bar: f64,
     zeta_bar: f64,
-    /// `|ζ̄₀| = ‖Âᵀb‖`, clamped positive; the reference for relative NE residuals.
+    /// `|ζ̄₀| = α₁ = ‖Âᵀu₁‖`, positive since a zero `α₁` never reaches the recurrence.
     zeta0: f64,
 }
 
@@ -165,15 +175,10 @@ impl NormalEqQr {
         }
     }
 
-    /// `|ζ̄ₖ|` — running estimate of `‖Aᵀ r_k‖` (Fong & Saunders).
+    /// `|ζ̄ₖ|` — running estimate of `‖Aᵀ r_k‖` for `u₁` (Fong & Saunders).
     fn normar(&self) -> f64 {
         self.zeta_bar.abs()
     }
-
-    /// `|ζ̄ₖ| / |ζ̄₀|` — normal-equation residual relative to `‖Aᵀb‖`; the `ζ̄₀` clamp guards it.
-    fn relative_normar(&self) -> f64 {
-        self.normar() / self.zeta0
-    }
 }
 
 /// Fong & Saunders' third rotation chain `Q̃` (§3.4, as in SciPy's `lsmr`), tracking `‖r_k‖`.
@@ -223,6 +228,8 @@ fn ratio(a: f64, b: f64) -> f64 {
 
 /// Vectors carried by the recurrence; `(h, h̄)` let `x` be built without the full `V_k` basis.
 pub(super) struct SolutionState {
+    /// `β₁`, restoring to `x` the scale the recurrence divided out of `b`.
+    rhs_norm: f64,
     x: Vec<f64>,
     h: Vec<f64>,
     h_bar: Vec<f64>,
@@ -230,8 +237,9 @@ pub(super) struct SolutionState {
 
 impl SolutionState {
     /// Initialize from the first normalized basis vector: `h₁ = v₁`, `x = 0`, `h̄₀ = 0`.
-    pub(super) fn init(v1: &[f64]) -> Self {
+    pub(super) fn init(v1: &[f64], rhs_norm: f64) -> Self {
         Self {
+            rhs_norm,
             x: vec![0.0; v1.len()],
             h: v1.to_vec(),
             h_bar: vec![0.0; v1.len()],
@@ -241,7 +249,7 @@ impl SolutionState {
     /// One `(x, h, h̄)` step; `v` must be normalized `v_{k+1}` and `prev` carries `(ρ, ρ̄)_{k-1}`.
     pub(super) fn update(&mut self, v: &[f64], curr: RotationStep, prev: RotationStep) {
         // One diagonal at a time: `ρρ̄` scales with `‖A‖²` and leaves the double range first.
-        let t_x = ratio(ratio(curr.zeta, curr.rho), curr.rho_bar);
+        let t_x = self.rhs_norm * ratio(ratio(curr.zeta, curr.rho), curr.rho_bar);
         let t_hbar = ratio(curr.theta_bar, prev.rho) * ratio(curr.rho, prev.rho_bar);
         let t_h = ratio(curr.theta_new, curr.rho);
 
@@ -308,7 +316,7 @@ impl ConvergenceCriteria {
     pub(super) fn start(self, alpha1: f64) -> ConvergenceState {
         ConvergenceState {
             criteria: self,
-            a_norm_sq: alpha1 * alpha1,
+            a_norm: alpha1,
         }
     }
 
@@ -323,7 +331,7 @@ impl ConvergenceCriteria {
     }
 
     /// An annihilated `α₁` certifies via the residual, or via a lower `‖A‖` bound that overstates.
-    pub(super) fn corroborates(&self, normr: f64, normar: f64, a_norm_below: f64) -> bool {
+    pub(super) fn corroborates(&self, normr: f64, normar: Magnitude, a_norm_below: f64) -> bool {
         normr <= self.residual_gap()
             || backward_error(normar, a_norm_below, normr) <= CERTIFICATION_SLACK * self.rel_tol
     }
@@ -332,27 +340,23 @@ impl ConvergenceCriteria {
 /// Mutable convergence observations for one LSMR run.
 pub(super) struct ConvergenceState {
     criteria: ConvergenceCriteria,
-    a_norm_sq: f64,
+    /// `‖A‖_F` estimate, never squared: `‖A‖²` leaves the double range from `‖A‖ ≈ 1e±154`.
+    a_norm: f64,
 }
 
 impl ConvergenceState {
-    /// Fold a fresh bidiagonal step into the `‖A‖_F²` estimate.
     pub(super) fn observe(&mut self, s: BidiagStep) {
-        self.a_norm_sq += s.alpha * s.alpha + s.beta * s.beta;
+        self.a_norm = self.a_norm.hypot(s.alpha).hypot(s.beta);
     }
 
     /// Check both stop criteria against the current scalar state.
     pub(super) fn check(&self, r: &LsmrRecurrenceState) -> Stop {
-        let residual = r.residual_estimate();
-        if residual <= self.criteria.abs_tol {
+        if r.residual_estimate() <= self.criteria.abs_tol {
             return Stop::ResidualTolerance;
         }
-        // The stream's own backward error, from its accumulated `‖A‖_F` estimate.
-        let ratio = backward_error(
-            r.normal_eq_residual_estimate(),
-            self.a_norm_sq.sqrt(),
-            residual,
-        );
+        // The stream's own backward error; `β₁` cancels, so the `u₁` form stays in range.
+        let normar = Magnitude::from(r.normal_eq_qr.normar());
+        let ratio = backward_error(normar, self.a_norm, r.unit_residual());
         if ratio <= self.criteria.rel_tol {
             return Stop::NormalEquationTolerance;
         }
@@ -361,17 +365,12 @@ impl ConvergenceState {
 }
 
 /// `‖Aᵀr‖ / (‖A‖‖r‖)`, refusing outright on a denominator carrying no information: clamping one
-/// up would flatter the ratio into certifying an unsolved stop. The product is the accurate form,
-/// rounding once; only where it leaves the float range does dividing in turn beat it.
-pub(super) fn backward_error(normar: f64, a_norm: f64, residual: f64) -> f64 {
+/// up would flatter the ratio into certifying an unsolved stop.
+pub(super) fn backward_error(normar: Magnitude, a_norm: f64, residual: f64) -> f64 {
     if !(a_norm > 0.0 && a_norm.is_finite() && residual > 0.0 && residual.is_finite()) {
         return f64::INFINITY;
     }
-    let denominator = a_norm * residual;
-    if denominator > 0.0 && denominator.is_finite() {
-        return normar / denominator;
-    }
-    normar / a_norm / residual
+    (normar / Magnitude::product(a_norm, residual)).to_f64()
 }
 
 /// Collapsed recurrences miss by orders of magnitude; the slack absorbs ordinary estimate drift.
diff --git a/crates/schwarz-precond/src/lsmr/tests/audit.rs b/crates/schwarz-precond/src/lsmr/tests/audit.rs
index 8bc96ba1..ec7b00c1 100644
--- a/crates/schwarz-precond/src/lsmr/tests/audit.rs
+++ b/crates/schwarz-precond/src/lsmr/tests/audit.rs
@@ -2,6 +2,7 @@
 use rstest::rstest;
 
 use super::super::bidiag::{BidiagStep, Bidiagonalization};
+use super::super::magnitude::Magnitude;
 use super::super::recurrence::ConvergenceCriteria;
 use super::super::{
     lsmr_from_bidiag, mlsmr, EscalationHandler, EscalationPolicy, LsmrResult, LsmrStopReason,
@@ -160,15 +161,15 @@ fn a_refused_stop_reports_its_iterates_normal_equation_residual(#[values(1.0, 4.
 #[rstest]
 #[case::usable(4.0, 0.25)]
 #[case::zero(0.0, 0.5)]
-#[case::negative(-1.0, 0.5)]
 #[case::overflowed(f64::INFINITY, 0.5)]
+#[case::nan(f64::NAN, 0.5)]
 fn a_report_only_moves_to_a_usable_reference(#[case] metric: f64, #[case] expected: f64) {
     let step1 = BidiagStep {
         alpha: 2.0,
         beta: 1.0,
     };
     assert_eq!(
-        NormalEqReference::warm(metric, step1).relative(1.0),
+        NormalEqReference::warm(Magnitude::from(metric), step1).relative(Magnitude::from(1.0)),
         expected
     );
 }
@@ -195,9 +196,11 @@ fn the_backward_error_survives_the_ends_of_the_range(
     #[case] residual: f64,
     #[case] expected: f64,
 ) {
-    let ratio = super::super::recurrence::backward_error(normar, a_norm, residual);
+    let ratio = super::super::recurrence::backward_error(Magnitude::from(normar), a_norm, residual);
     if expected.is_infinite() {
         assert_eq!(ratio, expected);
+    } else if expected < f64::MIN_POSITIVE {
+        assert_eq!(ratio.to_bits(), expected.to_bits());
     } else {
         assert!(
             (ratio - expected).abs() <= 1e-3 * expected.max(f64::MIN_POSITIVE),
diff --git a/crates/schwarz-precond/src/lsmr/tests/solve.rs b/crates/schwarz-precond/src/lsmr/tests/solve.rs
index 9e6c5099..82fa5fae 100644
--- a/crates/schwarz-precond/src/lsmr/tests/solve.rs
+++ b/crates/schwarz-precond/src/lsmr/tests/solve.rs
@@ -528,18 +528,6 @@ fn a_small_normed_design_still_gets_its_solution_update() {
     assert!((r.x[0] / 1e4 - 1.0).abs() < 1e-9, "{:?}", r.x);
 }
 
-/// `ρρ̄` scales with `‖A‖²`, so it leaves the double range long before `‖A‖` itself does.
-#[rstest]
-#[case::underflowing(1e-163)]
-#[case::overflowing(1e160)]
-fn an_extreme_scale_design_still_gets_its_solution_update(#[case] s: f64) {
-    let r = lsmr(&DiagOp(vec![s, 2.0 * s]), &[1.0, 1.0], 1e-10, 50, None).expect("extreme solve");
-    assert!(r.converged, "stop_reason: {:?}", r.stop_reason);
-    for (xi, want) in r.x.iter().zip([1.0 / s, 0.5 / s]) {
-        assert!((xi / want - 1.0).abs() < 1e-12, "{:?}", r.x);
-    }
-}
-
 /// A budget stop's `‖r_k‖` is the returned iterate's own, not LSQR's `|φ̄_k|`.
 #[rstest]
 fn the_residual_estimate_is_the_iterates_own(#[values(2, 4, 8)] maxiter: usize) {
@@ -559,3 +547,57 @@ fn the_residual_estimate_is_the_iterates_own(#[values(2, 4, 8)] maxiter: usize)
         r.residual_norm
     );
 }
+
+/// Power-of-two scaling is exact, so an extreme `‖A‖` or `‖b‖` must reproduce the unit-scale solve.
+#[rstest]
+#[case::small_design(-560, 0)]
+#[case::small_product(-560, -470)]
+#[case::large_design(530, 0)]
+#[case::large_product(530, 500)]
+fn a_power_of_two_scaling_reproduces_the_unit_scale_solve(
+    #[case] design_exp: i32,
+    #[case] rhs_exp: i32,
+    #[values(false)] metric: bool,
+) {
+    let solve = |op: &DenseOp, b: &[f64]| {
+        let identity = IdentityOp { n: op.cols };
+        match metric {
+            false => lsmr(op, b, 1e-4, 100, None),
+            true => mlsmr(op, b, &identity, 1e-4, 100, MlsmrOptions::default()),
+        }
+        .expect("scaled solve")
+    };
+    // Conditioned well enough that one rounding apart stays far below the checked agreement.
+    let op = DenseOp::vandermonde(30, 6);
+    // A sawtooth no low-degree polynomial fits keeps the residual far above the tolerance.
+    let b: Vec<f64> = (0..op.rows)
+        .map(|i| (1.0 + i as f64 / 29.0).ln() + 0.1 * (-1f64).powi(i as i32))
+        .collect();
+    let unit = solve(&op, &b);
+    assert_eq!(unit.stop_reason, LsmrStopReason::NormalEquationTolerance);
+
+    let (s, t) = (2f64.powi(design_exp), 2f64.powi(rhs_exp));
+    let scaled_op = DenseOp {
+        data: op.data.iter().map(|a| a * s).collect(),
+        ..op
+    };
+    let scaled_b: Vec<f64> = b.iter().map(|bi| bi * t).collect();
+    let r = solve(&scaled_op, &scaled_b);
+
+    assert_eq!(
+        (r.stop_reason, r.iterations),
+        (unit.stop_reason, unit.iterations)
+    );
+    assert!(r.converged);
+    assert!((r.residual_norm / t / unit.residual_norm - 1.0).abs() < 1e-10);
+    assert!(
+        (r.normal_eq_residual / unit.normal_eq_residual - 1.0).abs() < 1e-6,
+        "{:e} vs {:e}",
+        r.normal_eq_residual,
+        unit.normal_eq_residual
+    );
+    let x_max = unit.x.iter().fold(0.0f64, |m, x| m.max(x.abs()));
+    for (xi, ui) in r.x.iter().zip(&unit.x) {
+        assert!((xi * s / t - ui).abs() < 1e-8 * x_max, "{:?}", r.x);
+    }
+}

From 4c1ab6dafcfe6e93698cdf120d2d80357c61932a Mon Sep 17 00:00:00 2001
From: Kristof Schroeder <kristof_schroeder@web.de>
Date: Wed, 23 Sep 2026 19:17:55 +0200
Subject: [PATCH 2/4] =?UTF-8?q?refactor(lsmr):=20drop=20Magnitude's=20unus?=
 =?UTF-8?q?ed=20ordering,=20compute=20=E2=80=96r=E2=80=96=20once=20per=20c?=
 =?UTF-8?q?heck?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`NormalEqReference::relative` was the only user of `PartialOrd` and
`ZERO`; a zero-or-normal reference tests the same with `is_normal`.
`Mul` had one caller and is inlined into `product`. `check` evaluates
the unit residual once for both stop tests, and the scaling sweep drops
its never-run `mlsmr` arm. Bitwise identical on the probe set.
---
 crates/schwarz-precond/src/lsmr.rs            |  2 +-
 crates/schwarz-precond/src/lsmr/magnitude.rs  | 54 +++++--------------
 crates/schwarz-precond/src/lsmr/recurrence.rs |  8 +--
 .../schwarz-precond/src/lsmr/tests/solve.rs   | 10 +---
 4 files changed, 20 insertions(+), 54 deletions(-)

diff --git a/crates/schwarz-precond/src/lsmr.rs b/crates/schwarz-precond/src/lsmr.rs
index 532e961a..13933d4b 100644
--- a/crates/schwarz-precond/src/lsmr.rs
+++ b/crates/schwarz-precond/src/lsmr.rs
@@ -366,7 +366,7 @@ impl NormalEqReference {
 
     /// A vanished reference divides as the smallest normal instead.
     fn relative(self, estimate: Magnitude) -> f64 {
-        let reference = if self.0 > Magnitude::ZERO {
+        let reference = if self.0.is_normal() {
             self.0
         } else {
             Magnitude::from(f64::MIN_POSITIVE)
diff --git a/crates/schwarz-precond/src/lsmr/magnitude.rs b/crates/schwarz-precond/src/lsmr/magnitude.rs
index 0e393480..ce8b42e5 100644
--- a/crates/schwarz-precond/src/lsmr/magnitude.rs
+++ b/crates/schwarz-precond/src/lsmr/magnitude.rs
@@ -1,31 +1,28 @@
 //! Norms as `m · 2^e`, for comparing gradients of problems whose scales share no range.
 
-use std::cmp::Ordering;
-use std::ops::{Div, Mul};
+use std::ops::Div;
 
 /// A nonnegative `m · 2^e`, `m ∈ [1, 2)`, zero or non-finite, for expressions a few factors deep.
-#[derive(Clone, Copy, Debug, PartialEq)]
+#[derive(Clone, Copy)]
 pub(super) struct Magnitude {
     m: f64,
     e: i32,
 }
 
 impl Magnitude {
-    pub(super) const ZERO: Self = Self { m: 0.0, e: 0 };
-
     /// `a · b`, exact in exponent where the `f64` product would over- or underflow.
     pub(super) fn product(a: f64, b: f64) -> Self {
-        Self::from(a) * Self::from(b)
+        let (a, b) = (Self::from(a), Self::from(b));
+        Self::normalized(a.m * b.m, a.e + b.e)
     }
 
     /// The nearest `f64`, saturating to 0 or ∞; a subnormal result rounds a second time.
     pub(super) fn to_f64(self) -> f64 {
-        if !self.m.is_normal() || self.e >= f64::MAX_EXP {
-            return if self.m.is_normal() {
-                f64::INFINITY
-            } else {
-                self.m
-            };
+        if !self.m.is_normal() {
+            return self.m;
+        }
+        if self.e >= f64::MAX_EXP {
+            return f64::INFINITY;
         }
         if self.e >= f64::MIN_EXP - 1 {
             return self.m * pow2(self.e);
@@ -67,13 +64,6 @@ impl From<f64> for Magnitude {
     }
 }
 
-impl Mul for Magnitude {
-    type Output = Self;
-    fn mul(self, rhs: Self) -> Self {
-        Self::normalized(self.m * rhs.m, self.e + rhs.e)
-    }
-}
-
 impl Div for Magnitude {
     type Output = Self;
     fn div(self, rhs: Self) -> Self {
@@ -81,16 +71,6 @@ impl Div for Magnitude {
     }
 }
 
-impl PartialOrd for Magnitude {
-    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
-        if self.m.is_normal() && other.m.is_normal() {
-            return Some(self.e.cmp(&other.e).then(self.m.total_cmp(&other.m)));
-        }
-        // 0, ∞ and NaN order against any `m ∈ [1, 2)` as they do against its value.
-        self.m.partial_cmp(&other.m)
-    }
-}
-
 /// `2^e` for a normal exponent, built from its bits since `powi` underflows on the way there.
 fn pow2(e: i32) -> f64 {
     debug_assert!((f64::MIN_EXP - 1..f64::MAX_EXP).contains(&e));
@@ -144,16 +124,10 @@ mod tests {
     }
 
     #[test]
-    fn the_order_is_the_value_order() {
-        let tiny = Magnitude::product(1e-300, 1e-300);
-        let huge = Magnitude::product(1e300, 1e300);
-        assert!(Magnitude::ZERO < tiny && tiny < Magnitude::from(1.0));
-        assert!(Magnitude::from(1.0) < huge && huge < Magnitude::from(f64::INFINITY));
-        assert!(Magnitude::from(1.5) > Magnitude::from(1.25));
-        assert_eq!(tiny.to_f64(), 0.0);
-        assert_eq!(huge.to_f64(), f64::INFINITY);
-        assert!(Magnitude::from(f64::NAN)
-            .partial_cmp(&Magnitude::ZERO)
-            .is_none());
+    fn an_out_of_range_value_saturates() {
+        assert!(Magnitude::product(1e-300, 1e-300).is_normal());
+        assert_eq!(Magnitude::product(1e-300, 1e-300).to_f64(), 0.0);
+        assert_eq!(Magnitude::product(1e300, 1e300).to_f64(), f64::INFINITY);
+        assert!(Magnitude::from(f64::NAN).to_f64().is_nan());
     }
 }
diff --git a/crates/schwarz-precond/src/lsmr/recurrence.rs b/crates/schwarz-precond/src/lsmr/recurrence.rs
index 150ffec7..b80b2104 100644
--- a/crates/schwarz-precond/src/lsmr/recurrence.rs
+++ b/crates/schwarz-precond/src/lsmr/recurrence.rs
@@ -351,12 +351,13 @@ impl ConvergenceState {
 
     /// Check both stop criteria against the current scalar state.
     pub(super) fn check(&self, r: &LsmrRecurrenceState) -> Stop {
-        if r.residual_estimate() <= self.criteria.abs_tol {
+        let unit_residual = r.unit_residual();
+        if r.rhs_norm * unit_residual <= self.criteria.abs_tol {
             return Stop::ResidualTolerance;
         }
         // The stream's own backward error; `β₁` cancels, so the `u₁` form stays in range.
         let normar = Magnitude::from(r.normal_eq_qr.normar());
-        let ratio = backward_error(normar, self.a_norm, r.unit_residual());
+        let ratio = backward_error(normar, self.a_norm, unit_residual);
         if ratio <= self.criteria.rel_tol {
             return Stop::NormalEquationTolerance;
         }
@@ -364,8 +365,7 @@ impl ConvergenceState {
     }
 }
 
-/// `‖Aᵀr‖ / (‖A‖‖r‖)`, refusing outright on a denominator carrying no information: clamping one
-/// up would flatter the ratio into certifying an unsolved stop.
+/// `‖Aᵀr‖ / (‖A‖‖r‖)`; an uninformative denominator refuses, since clamping it would certify.
 pub(super) fn backward_error(normar: Magnitude, a_norm: f64, residual: f64) -> f64 {
     if !(a_norm > 0.0 && a_norm.is_finite() && residual > 0.0 && residual.is_finite()) {
         return f64::INFINITY;
diff --git a/crates/schwarz-precond/src/lsmr/tests/solve.rs b/crates/schwarz-precond/src/lsmr/tests/solve.rs
index 82fa5fae..1c25cbb4 100644
--- a/crates/schwarz-precond/src/lsmr/tests/solve.rs
+++ b/crates/schwarz-precond/src/lsmr/tests/solve.rs
@@ -557,16 +557,8 @@ fn the_residual_estimate_is_the_iterates_own(#[values(2, 4, 8)] maxiter: usize)
 fn a_power_of_two_scaling_reproduces_the_unit_scale_solve(
     #[case] design_exp: i32,
     #[case] rhs_exp: i32,
-    #[values(false)] metric: bool,
 ) {
-    let solve = |op: &DenseOp, b: &[f64]| {
-        let identity = IdentityOp { n: op.cols };
-        match metric {
-            false => lsmr(op, b, 1e-4, 100, None),
-            true => mlsmr(op, b, &identity, 1e-4, 100, MlsmrOptions::default()),
-        }
-        .expect("scaled solve")
-    };
+    let solve = |op: &DenseOp, b: &[f64]| lsmr(op, b, 1e-4, 100, None).expect("scaled solve");
     // Conditioned well enough that one rounding apart stays far below the checked agreement.
     let op = DenseOp::vandermonde(30, 6);
     // A sawtooth no low-degree polynomial fits keeps the residual far above the tolerance.

From 4bd25ce8e5b4d5088f79e68bd91908e7e979687a Mon Sep 17 00:00:00 2001
From: Kristof Schroeder <kristof_schroeder@web.de>
Date: Thu, 24 Sep 2026 14:09:50 +0200
Subject: [PATCH 3/4] =?UTF-8?q?fix(lsmr):=20form=20the=20solution=20step's?=
 =?UTF-8?q?=20=CE=B2=E2=82=81=CE=B6/(=CF=81=CF=81=CC=84)=20in=20exact=20ex?=
 =?UTF-8?q?ponents?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`ζ/(ρρ̄)` on the unit RHS can underflow to zero where `β₁ζ/(ρρ̄)` is a normal
double, so `x` stayed 0 while the recurrence reported a converged stop.

Co-Authored-By: Claude <noreply@anthropic.com>
---
 crates/schwarz-precond/src/lsmr/recurrence.rs | 10 ++++++--
 .../schwarz-precond/src/lsmr/tests/solve.rs   | 24 +++++++++++++++++++
 2 files changed, 32 insertions(+), 2 deletions(-)

diff --git a/crates/schwarz-precond/src/lsmr/recurrence.rs b/crates/schwarz-precond/src/lsmr/recurrence.rs
index b80b2104..9c5437ae 100644
--- a/crates/schwarz-precond/src/lsmr/recurrence.rs
+++ b/crates/schwarz-precond/src/lsmr/recurrence.rs
@@ -248,8 +248,14 @@ impl SolutionState {
 
     /// One `(x, h, h̄)` step; `v` must be normalized `v_{k+1}` and `prev` carries `(ρ, ρ̄)_{k-1}`.
     pub(super) fn update(&mut self, v: &[f64], curr: RotationStep, prev: RotationStep) {
-        // One diagonal at a time: `ρρ̄` scales with `‖A‖²` and leaves the double range first.
-        let t_x = self.rhs_norm * ratio(ratio(curr.zeta, curr.rho), curr.rho_bar);
+        // `ρρ̄` scales with `‖A‖²`, and `ζ/(ρρ̄)` can underflow where `β₁ζ/(ρρ̄)` does not.
+        let t_x = if curr.rho == 0.0 || curr.rho_bar == 0.0 {
+            0.0
+        } else {
+            let unsigned = Magnitude::product(self.rhs_norm, curr.zeta.abs())
+                / Magnitude::product(curr.rho, curr.rho_bar);
+            curr.zeta.signum() * unsigned.to_f64()
+        };
         let t_hbar = ratio(curr.theta_bar, prev.rho) * ratio(curr.rho, prev.rho_bar);
         let t_h = ratio(curr.theta_new, curr.rho);
 
diff --git a/crates/schwarz-precond/src/lsmr/tests/solve.rs b/crates/schwarz-precond/src/lsmr/tests/solve.rs
index 1c25cbb4..786d1ed4 100644
--- a/crates/schwarz-precond/src/lsmr/tests/solve.rs
+++ b/crates/schwarz-precond/src/lsmr/tests/solve.rs
@@ -528,6 +528,30 @@ fn a_small_normed_design_still_gets_its_solution_update() {
     assert!((r.x[0] / 1e4 - 1.0).abs() < 1e-9, "{:?}", r.x);
 }
 
+/// `ζ/(ρρ̄)` rounds to zero below `2^-1074` although `β₁` scales it back to a normal `x`.
+#[rstest]
+#[case::overflowing_rhs(960, 1023)]
+#[case::dominant_miss(23, 80)]
+fn the_solution_update_survives_its_unit_rhs_underflowing(
+    #[case] fit_exp: i32,
+    #[case] miss_exp: i32,
+) {
+    let op = DenseOp {
+        rows: 2,
+        cols: 1,
+        data: vec![2f64.powi(1023), 0.0],
+    };
+    let b = [2f64.powi(fit_exp), 2f64.powi(miss_exp)];
+    let r = lsmr(&op, &b, 1e-10, 10, None).expect("extreme solve");
+    assert!(r.converged, "stop_reason: {:?}", r.stop_reason);
+    let want = 2f64.powi(fit_exp - 1023);
+    assert!(
+        (r.x[0] / want - 1.0).abs() < 1e-12,
+        "{:e} vs {want:e}",
+        r.x[0]
+    );
+}
+
 /// A budget stop's `‖r_k‖` is the returned iterate's own, not LSQR's `|φ̄_k|`.
 #[rstest]
 fn the_residual_estimate_is_the_iterates_own(#[values(2, 4, 8)] maxiter: usize) {

From 13e937aa1bd824720f7d67f180288c213298fd90 Mon Sep 17 00:00:00 2001
From: Kristof Schroeder <kristof_schroeder@web.de>
Date: Thu, 24 Sep 2026 14:10:09 +0200
Subject: [PATCH 4/4] refactor(lsmr): drop NormalEqReference's unreachable
 zero-reference fallback
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`Magnitude::product` of two positive finite norms never vanishes, and a zero α₁
returns before any report; also corrects the `warm` and `pow2` docs.

Co-Authored-By: Claude <noreply@anthropic.com>
---
 crates/schwarz-precond/src/lsmr.rs             | 11 +++--------
 crates/schwarz-precond/src/lsmr/magnitude.rs   |  2 +-
 crates/schwarz-precond/src/lsmr/tests/solve.rs |  2 +-
 3 files changed, 5 insertions(+), 10 deletions(-)

diff --git a/crates/schwarz-precond/src/lsmr.rs b/crates/schwarz-precond/src/lsmr.rs
index 13933d4b..6d8ebad8 100644
--- a/crates/schwarz-precond/src/lsmr.rs
+++ b/crates/schwarz-precond/src/lsmr.rs
@@ -355,7 +355,7 @@ impl NormalEqReference {
         Self(Magnitude::product(step1.alpha, step1.beta))
     }
 
-    /// A reference that carries no information leaves the stream's own `ζ̄₀` to divide.
+    /// A reference that carries no information falls back to the cold `α₁β₁`.
     fn warm(metric: Magnitude, step1: BidiagStep) -> Self {
         if metric.is_normal() {
             Self(metric)
@@ -364,14 +364,9 @@ impl NormalEqReference {
         }
     }
 
-    /// A vanished reference divides as the smallest normal instead.
     fn relative(self, estimate: Magnitude) -> f64 {
-        let reference = if self.0.is_normal() {
-            self.0
-        } else {
-            Magnitude::from(f64::MIN_POSITIVE)
-        };
-        (estimate / reference).to_f64()
+        debug_assert!(self.0.is_normal(), "a zero α₁ returns before any report");
+        (estimate / self.0).to_f64()
     }
 }
 
diff --git a/crates/schwarz-precond/src/lsmr/magnitude.rs b/crates/schwarz-precond/src/lsmr/magnitude.rs
index ce8b42e5..7c2b7204 100644
--- a/crates/schwarz-precond/src/lsmr/magnitude.rs
+++ b/crates/schwarz-precond/src/lsmr/magnitude.rs
@@ -71,7 +71,7 @@ impl Div for Magnitude {
     }
 }
 
-/// `2^e` for a normal exponent, built from its bits since `powi` underflows on the way there.
+/// `2^e` for a normal exponent, built from its bits since `powi`'s precision is unspecified.
 fn pow2(e: i32) -> f64 {
     debug_assert!((f64::MIN_EXP - 1..f64::MAX_EXP).contains(&e));
     f64::from_bits(((e + 1023) as u64) << 52)
diff --git a/crates/schwarz-precond/src/lsmr/tests/solve.rs b/crates/schwarz-precond/src/lsmr/tests/solve.rs
index 786d1ed4..7961e555 100644
--- a/crates/schwarz-precond/src/lsmr/tests/solve.rs
+++ b/crates/schwarz-precond/src/lsmr/tests/solve.rs
@@ -587,7 +587,7 @@ fn a_power_of_two_scaling_reproduces_the_unit_scale_solve(
     let op = DenseOp::vandermonde(30, 6);
     // A sawtooth no low-degree polynomial fits keeps the residual far above the tolerance.
     let b: Vec<f64> = (0..op.rows)
-        .map(|i| (1.0 + i as f64 / 29.0).ln() + 0.1 * (-1f64).powi(i as i32))
+        .map(|i| (1.0 + i as f64 / (op.rows - 1) as f64).ln() + 0.1 * (-1f64).powi(i as i32))
         .collect();
     let unit = solve(&op, &b);
     assert_eq!(unit.stop_reason, LsmrStopReason::NormalEquationTolerance);