docs(#327): link the fixed-thread-count reproducibility contract, and stop asserting cross-count bitwise equality

faer/src/lib.rs

@@ -180,6 +180,15 @@
 //! The operator overloads (`&a * &b`, `&a + &b`, ...) have no room for the argument and always run
 //! sequentially; use [`linalg::matmul::matmul`] for a parallel product.
 //!
+//! Results are bitwise reproducible on the same machine and build **at a fixed thread count**,
+//! and only there: a different `Par` width is not promised to agree bit for bit. The width can
+//! change both the path a call takes (the cost model gates the sequential and the parallel paths
+//! on the problem and the width it is given) and how the work is partitioned (some routines split
+//! by the participant count itself, such as the trailing update of the tridiagonalization), so two
+//! widths can group the same sum differently. See
+//! [docs/reproducibility.md](https://github.com/tensor4all/faer/blob/main/docs/reproducibility.md)
+//! (tensor4all/faer#327).
+//!
 //! with the `threadteam` feature (on by default) a team is a reservation the caller owns: its members spin while the reservation is held,
 //! instead of being woken for every region, and the token hands each region the authority to
 //! place work on them. That is what makes a sequence of short regions cheap: the

faer/src/linalg/matmul/mod.rs

@@ -2456,6 +2456,11 @@ mod tests {
 		check_uninit_native(|x, y| c64::new(x, y), |a: c64, b| (a - b).norm(), 1e-14);
 	}
 
+	/// The uninitialized entry matches the replace entry on a team: the reference is built with
+	/// the same team and compared with `==`, which treats `+0.0` and `-0.0` as equal, so this is
+	/// representation equality at one width rather than a bitwise claim (bitwise agreement is
+	/// promised at a fixed thread count only, `docs/reproducibility.md`, tensor4all/faer#327).
+	/// The serial product is compared within a relative tolerance.
 	#[cfg(feature = "threadteam")]
 	#[test]
 	fn test_matmul_with_conj_uninit_team() {
@@ -2464,18 +2469,51 @@ mod tests {
 				let lhs = Mat::<f64>::from_fn(130, 70, |i, j| ((i * 7 + j * 3) % 11) as f64 - 5.0);
 				let rhs = Mat::<f64>::from_fn(70, 90, |i, j| ((i * 5 + j) % 13) as f64 - 6.0);
 				let mut expected = Mat::<f64>::zeros(130, 90);
-				matmul_with_conj(expected.as_mut(), Accum::Replace, lhs.as_ref(), Conj::No, rhs.as_ref(), Conj::No, 1.0, Par::Seq);
+				matmul_with_conj(expected.as_mut(), Accum::Replace, lhs.as_ref(), Conj::No, rhs.as_ref(), Conj::No, 1.0, par.reborrow());
 				let mut data = alloc::vec![MaybeUninit::new(f64::NAN); 130 * 90];
 				let dst = MatUninitMut::from_column_major_slice_mut(&mut data, 130, 90);
 				let result = matmul_with_conj_uninit(dst, lhs.as_ref(), Conj::No, rhs.as_ref(), Conj::No, 1.0, par.reborrow());
 				assert!(result == expected);
+				let mut sequential = Mat::<f64>::zeros(130, 90);
+				matmul_with_conj(sequential.as_mut(), Accum::Replace, lhs.as_ref(), Conj::No, rhs.as_ref(), Conj::No, 1.0, Par::Seq);
+				// Scaled by the product itself, so the bound states a relative error: a
+				// near-zero entry is compared against the scale of the whole matrix rather
+				// than against a floor that would admit a much larger relative error. An
+				// all-zero reference has no scale to be relative to; it is outside this test,
+				// and the assert below rejects it before any comparison.
+				let bound_of = |m: MatRef<'_, f64>| {
+					let scale = (0..90)
+						.flat_map(|j| (0..130).map(move |i| (i, j)))
+						.fold(0.0f64, |scale, (i, j)| scale.max(m[(i, j)].abs()));
+					assert!(scale > 0.0, "a reference of zeros is outside this test");
+					1e-12 * scale
+				};
+				let bound = bound_of(sequential.as_ref());
+				for j in 0..90 {
+					for i in 0..130 {
+						let diff = (result[(i, j)] - sequential[(i, j)]).abs();
+						::core::assert!(diff <= bound, "{i}, {j}: {diff} against {bound}");
+					}
+				}
+				// The same product shrunk by 1e-20, where the largest entry is 1.01e-18:
+				// an all-zero result must exceed the bound, which it cannot do under a
+				// `max(1, ...)` floor (that bound is 1e-12 whatever the entries are).
+				let tiny = Mat::<f64>::from_fn(130, 90, |i, j| 1e-20 * sequential[(i, j)]);
+				let tiny_bound = bound_of(tiny.as_ref());
+				let zeros = (0..90)
+					.flat_map(|j| (0..130).map(move |i| (i, j)))
+					.fold(0.0f64, |d, (i, j)| d.max(tiny[(i, j)].abs()));
+				::core::assert!(
+					zeros > tiny_bound,
+					"{zeros} against {tiny_bound}: the bound admits an all-zero product"
+				);
 			}
 			{
 				let from = |x: f64| <fx128 as ComplexFieldExt>::from_f64(x);
 				let lhs = Mat::<fx128>::from_fn(9, 6, |i, j| from((i + 2 * j) as f64));
 				let rhs = Mat::<fx128>::from_fn(6, 11, |i, j| from((3 * i + j) as f64));
 				let mut expected = Mat::<fx128>::zeros(9, 11);
-				matmul_with_conj(expected.as_mut(), Accum::Replace, lhs.as_ref(), Conj::No, rhs.as_ref(), Conj::No, from(1.0), Par::Seq);
+				matmul_with_conj(expected.as_mut(), Accum::Replace, lhs.as_ref(), Conj::No, rhs.as_ref(), Conj::No, from(1.0), par.reborrow());
 				let mut data = alloc::vec![MaybeUninit::<fx128>::uninit(); 9 * 11];
 				let dst = MatUninitMut::from_column_major_slice_mut(&mut data, 9, 11);
 				let result = matmul_with_conj_uninit(dst, lhs.as_ref(), Conj::No, rhs.as_ref(), Conj::No, from(1.0), par.reborrow());

faer/src/utils/team.rs

@@ -859,34 +859,131 @@ mod tests {
 		})
 	}
 
+	/// Folds one difference into the running maximum, refusing a non-finite one: `f64::max`
+	/// returns its other operand when one side is NaN, so an all-NaN output would satisfy any
+	/// tolerance.
+	fn max_finite(d: f64, x: f64) -> f64 {
+		assert!(x.is_finite(), "non-finite difference {x}");
+		d.max(x)
+	}
+
 	fn max_diff(a: MatRef<'_, f64>, b: MatRef<'_, f64>) -> f64 {
 		let mut d: f64 = 0.0;
 		for j in 0..a.ncols() {
 			for i in 0..a.nrows() {
-				d = d.max((a[(i, j)] - b[(i, j)]).abs());
+				let (av, bv) = (a[(i, j)], b[(i, j)]);
+				let x = (av - bv).abs();
+				// A non-finite difference is accepted only between two equal infinities: `q`
+				// below holds `+inf` at (3, 31) on both routes. A NaN never passes, not even
+				// two NaNs of the same bits, whose difference is itself NaN.
+				assert!(x.is_finite() || (av == bv && av.is_infinite()), "{i},{j}: {av} and {bv} differ by {x}");
+				// `x` is finite here, or both entries are the same infinity and `x` is NaN,
+				// which `f64::max` drops.
+				d = d.max(x);
 			}
 		}
 		d
 	}
 
+	/// The exception in `max_diff` is for the `+inf` of the 32x32 QR fixture, not for NaNs: two
+	/// NaNs of the same bits have a NaN difference, which `f64::max` would drop.
+	#[test]
+	#[should_panic(expected = "differ by NaN")]
+	fn max_diff_rejects_a_nan_on_both_sides() {
+		let nan = Mat::<f64>::full(1, 1, f64::NAN);
+		max_diff(nan.as_ref(), nan.as_ref());
+	}
+
+	#[test]
+	fn max_diff_accepts_two_equal_infinities() {
+		let inf = Mat::<f64>::full(1, 1, f64::INFINITY);
+		assert_eq!(max_diff(inf.as_ref(), inf.as_ref()), 0.0);
+	}
+
+	/// Relative comparison for the families below: `max_diff(a, b)` within `1e-10` of the largest
+	/// finite magnitude of the reference `a` (at least 1). Non-finite entries stay out of the
+	/// scale: the `+inf` of the 32x32 QR fixture would make the bound infinite and let any finite
+	/// error through. `max_diff` still requires such entries to be the same infinity on both sides.
+	fn assert_within_tolerance(a: MatRef<'_, f64>, b: MatRef<'_, f64>, what: &str) {
+		let scale = (0..a.ncols()).flat_map(|j| (0..a.nrows()).map(move |i| (i, j))).fold(1.0f64, |m, (i, j)| {
+			let v = a[(i, j)].abs();
+			if v.is_finite() { m.max(v) } else { m }
+		});
+		let d = max_diff(a, b);
+		assert!(d <= 1e-10 * scale, "{what}: {d} against {scale}");
+	}
+
+	/// The 32x32 no-pivoting QR used by `small_qr_on_a_team_is_the_sequential_one`; its `q` holds
+	/// `+inf` at (3, 31).
+	fn small_qr(mut par: Par<'_>) -> (Mat<f64>, Mat<f64>) {
+		use crate::dyn_stack::{MemBuffer, MemStack};
+		use crate::linalg::qr::no_pivoting::factor as qrf;
+		let n = 32;
+		let bs = qrf::recommended_block_size::<f64>(n, n);
+		let mut a = rand_mat(n, n, 11);
+		let mut q = Mat::<f64>::zeros(bs, n);
+		let mut mem = MemBuffer::new(qrf::qr_in_place_scratch::<f64>(n, n, bs, par.shape(), Default::default()));
+		qrf::qr_in_place(a.as_mut(), q.as_mut(), par.reborrow(), MemStack::new(&mut mem), Default::default());
+		(a, q)
+	}
+
+	/// The 200x200 column pivoting QR of the families test: the factor `a` and the Householder
+	/// coefficients `h`, which hold an infinity sentinel.
+	fn col_piv_qr(par: Par<'_>) -> Vec<Mat<f64>> {
+		use crate::dyn_stack::{MemBuffer, MemStack};
+		use crate::linalg;
+		let n = 200;
+		let mut a = rand_mat(n, n, 13);
+		let bs = linalg::qr::no_pivoting::factor::recommended_block_size::<f64>(n, n);
+		let mut h = Mat::<f64>::zeros(bs, n);
+		let (mut q, mut qi) = (vec![0usize; n], vec![0usize; n]);
+		let req = linalg::qr::col_pivoting::factor::qr_in_place_scratch::<usize, f64>(n, n, bs, par.shape(), Default::default());
+		linalg::qr::col_pivoting::factor::qr_in_place(a.as_mut(), h.as_mut(), &mut q, &mut qi, par, MemStack::new(&mut MemBuffer::new(req)), Default::default());
+		vec![a, h]
+	}
+
+	/// The comparison the families test applies, on its own column pivoting QR: `h` holds an
+	/// infinity, and a copy whose (0, 0) is off by `1e6` must still be rejected.
+	#[test]
+	#[should_panic(expected = "against")]
+	fn tolerance_rejects_a_corrupted_col_pivoting_factor() {
+		let h = col_piv_qr(Par::Seq).pop().unwrap();
+		let infinite = (0..h.ncols()).flat_map(|j| (0..h.nrows()).map(move |i| (i, j))).any(|(i, j)| h[(i, j)].is_infinite());
+		assert!(infinite, "the fixture lost its sentinel");
+		assert!(h[(0, 0)].is_finite());
+		let mut off = h.clone();
+		off[(0, 0)] += 1e6;
+		assert_within_tolerance(h.as_ref(), off.as_ref(), "column pivoting h with (0, 0) off by 1e6");
+	}
+
+	/// An infinity in the reference must not widen the tolerance: the QR fixture's `q` against a
+	/// copy whose (0, 0) is off by `1e6` is rejected although (3, 31) is `+inf` on both sides.
+	#[test]
+	#[should_panic(expected = "against")]
+	fn tolerance_ignores_an_infinite_reference_entry() {
+		let (_, q) = small_qr(Par::Seq);
+		assert_eq!(q[(3, 31)], f64::INFINITY, "the fixture lost its sentinel");
+		let mut off = q.clone();
+		off[(0, 0)] += 1e6;
+		assert_within_tolerance(q.as_ref(), off.as_ref(), "q with (0, 0) off by 1e6");
+	}
+
+	/// The guard the manual comparisons rely on, which no `max_diff` test reaches: a NaN must not
+	/// fold into the running maximum.
+	#[test]
+	#[should_panic(expected = "non-finite difference")]
+	fn max_finite_rejects_a_nan() {
+		max_finite(0.0, f64::NAN);
+	}
+
 	/// A factorization too small to pay for a section asks the gates of a team token as the
 	/// sequential path does: the same factors to the last bit, and no section opened. The block
 	/// apply decides once at its top (`householder.rs`, the floor of the whole apply) and runs its
 	/// products sequentially, instead of asking each product and each triangular step.
 	#[test]
 	fn small_qr_on_a_team_is_the_sequential_one() {
-		use crate::dyn_stack::{MemBuffer, MemStack};
-		use crate::linalg::qr::no_pivoting::factor as qrf;
 		let n = 32;
-		let a0 = rand_mat(n, n, 11);
-		let run = |mut par: Par<'_>| {
-			let bs = qrf::recommended_block_size::<f64>(n, n);
-			let mut a = a0.clone();
-			let mut q = Mat::<f64>::zeros(bs, n);
-			let mut mem = MemBuffer::new(qrf::qr_in_place_scratch::<f64>(n, n, bs, par.shape(), Default::default()));
-			qrf::qr_in_place(a.as_mut(), q.as_mut(), par.reborrow(), MemStack::new(&mut mem), Default::default());
-			(a, q)
-		};
+		let run = small_qr;
 		let (a_seq, q_seq) = run(Par::Seq);
 		let (a_team, q_team) = with_par(4, |par| {
 			stats::take();
@@ -899,6 +996,11 @@ mod tests {
 			assert_eq!(s.gemm_tiled + s.gemm_team_kernel, 0, "{s:?}");
 			r
 		});
+		// An exactly-zero difference is legitimate here because the team opened no section
+		// (`sections` and the GEMM counters above are zero), so both sides run the same
+		// one-participant partition. The `max_diff` form does not distinguish `+0.0` from
+		// `-0.0`; across thread counts only a tolerance is promised
+		// (`docs/reproducibility.md`).
 		assert_eq!(max_diff(a_seq.as_ref(), a_team.as_ref()), 0.0);
 		assert_eq!(max_diff(q_seq.as_ref(), q_team.as_ref()), 0.0);
 	}
@@ -1047,11 +1149,6 @@ mod tests {
 		let ns = StdDuration::from_nanos;
 		// cheap up to 4, prohibitive from 5 on
 		let row = [(2, ns(10)), (4, ns(10)), (5, StdDuration::from_secs(1))];
-		let tol = |a: MatRef<'_, f64>, b: MatRef<'_, f64>, what: &str| {
-			let scale = (0..a.ncols()).flat_map(|j| (0..a.nrows()).map(move |i| (i, j))).fold(1.0f64, |m, (i, j)| m.max(a[(i, j)].abs()));
-			let d = max_diff(a, b);
-			assert!(d <= 1e-10 * scale, "{what}: {d} against {scale}");
-		};
 		type Run = fn(Par<'_>) -> Vec<Mat<f64>>;
 		let families: [(&str, Run); 14] = [
 			("full pivoting lu", |par| {
@@ -1070,16 +1167,7 @@ mod tests {
 				linalg::lu::partial_pivoting::factor::lu_in_place(a.as_mut(), &mut p, &mut pi, par, MemStack::new(&mut MemBuffer::new(req)), Default::default());
 				vec![a]
 			}),
-			("column pivoting qr", |par| {
-				let n = 200;
-				let mut a = rand_mat(n, n, 13);
-				let bs = linalg::qr::no_pivoting::factor::recommended_block_size::<f64>(n, n);
-				let mut h = Mat::<f64>::zeros(bs, n);
-				let (mut q, mut qi) = (vec![0usize; n], vec![0usize; n]);
-				let req = linalg::qr::col_pivoting::factor::qr_in_place_scratch::<usize, f64>(n, n, bs, par.shape(), Default::default());
-				linalg::qr::col_pivoting::factor::qr_in_place(a.as_mut(), h.as_mut(), &mut q, &mut qi, par, MemStack::new(&mut MemBuffer::new(req)), Default::default());
-				vec![a, h]
-			}),
+			("column pivoting qr", col_piv_qr),
 			("hessenberg", |par| {
 				let n = 160;
 				let mut a = rand_mat(n, n, 14);
@@ -1207,7 +1295,7 @@ mod tests {
 			assert!(s.split[4] + s.run[4] > 0, "{name}: no split or section of 4 ({s:?})");
 			assert!(s.split.iter().skip(5).all(|&c| c == 0), "{name}: a split wider than 4 ({s:?})");
 			for (k, (a, b)) in want.iter().zip(got.iter()).enumerate() {
-				tol(a.as_ref(), b.as_ref(), &alloc::format!("{name} output {k}"));
+				assert_within_tolerance(a.as_ref(), b.as_ref(), &alloc::format!("{name} output {k}"));
 			}
 		}
 	}
@@ -1415,7 +1503,7 @@ mod tests {
 			let mut d: f64 = 0.0;
 			for j in 0..n {
 				for i in 0..m {
-					d = d.max((want[(i, j)] - got[(i, j)]).norm_sqr().sqrt());
+					d = max_finite(d, (want[(i, j)] - got[(i, j)]).norm_sqr().sqrt());
 				}
 			}
 			assert!(d < 1e-9);
@@ -1690,7 +1778,7 @@ mod tests {
 				let mut d: f64 = 0.0;
 				for j in 0..n {
 					for i in j..n {
-						d = d.max((want[(i, j)] - got[(i, j)]).abs());
+						d = max_finite(d, (want[(i, j)] - got[(i, j)]).abs());
 					}
 				}
 				assert!(d < 1e-8, "cholesky n = {n}: {d}");
@@ -1857,6 +1945,8 @@ mod tests {
 				matmul(got.as_mut(), Accum::Replace, a.as_ref(), b.as_ref(), 1.0, par);
 				let st = stats::take();
 				stats::enable(false);
+				// A team of one participant is the sequential partition, so an exactly-zero
+				// difference here is not a cross-thread-count claim (`docs/reproducibility.md`).
 				assert_eq!(max_diff(want.as_ref(), got.as_ref()), 0.0);
 				assert_eq!(st.gemm_team_kernel + st.gemm_tiled, 0, "{st:?}");
 			})
@@ -1889,7 +1979,7 @@ mod tests {
 			for j in 0..n {
 				for i in 0..n {
 					if (lo && i >= j) || (!lo && i <= j) {
-						d = d.max((m[(i, j)] - w[(i, j)]).abs());
+						d = max_finite(d, (m[(i, j)] - w[(i, j)]).abs());
 					}
 				}
 			}
@@ -1952,7 +2042,7 @@ mod tests {
 		let mut d: f64 = 0.0;
 		for j in 0..n {
 			for i in 0..n {
-				d = d.max((want[(i, j)] - got[(i, j)]).norm());
+				d = max_finite(d, (want[(i, j)] - got[(i, j)]).norm());
 			}
 		}
 		assert!(d < 1e-9, "{d}");
@@ -1981,7 +2071,7 @@ mod tests {
 		let mut d: f64 = 0.0;
 		for j in 0..n {
 			for i in j..n {
-				d = d.max((want[(i, j)] - got[(i, j)]).abs());
+				d = max_finite(d, (want[(i, j)] - got[(i, j)]).abs());
 			}
 		}
 		assert!(d < 1e-9, "{d}");
@@ -2038,7 +2128,7 @@ mod tests {
 							acc, alpha, conj, conj,
 						);
 					}
-					let d = want.iter().zip(&got).map(|(&w, &g)| err(w, g)).fold(0.0, f64::max);
+					let d = want.iter().zip(&got).map(|(&w, &g)| err(w, g)).fold(0.0, max_finite);
 					assert!(d < 1e-3, "read_dst {read_dst} conj {conj}: {d}");
 				}
 			}
@@ -2163,7 +2253,7 @@ mod tests {
 				let mut d = 0.0f64;
 				for j in 0..n {
 					for i in 0..m {
-						d = d.max(err(want[(i, j)], got[(i, j)]));
+						d = max_finite(d, err(want[(i, j)], got[(i, j)]));
 					}
 				}
 				assert!(d < 1e-3, "{beta:?}: {d}");
@@ -2196,7 +2286,7 @@ mod tests {
 			for j in 0..n {
 				for i in 0..n {
 					if (lower && i >= j) || (!lower && i <= j) {
-						d = d.max((x[(i, j)] - y[(i, j)]).abs());
+						d = max_finite(d, (x[(i, j)] - y[(i, j)]).abs());
 					}
 				}
 			}
@@ -2258,7 +2348,7 @@ mod tests {
 			let mut d: f64 = 0.0;
 			for j in 0..n {
 				for i in j..n {
-					d = d.max((x[(i, j)] - y[(i, j)]).abs());
+					d = max_finite(d, (x[(i, j)] - y[(i, j)]).abs());
 				}
 			}
 			d

faer/tests/common/mod.rs

@@ -22,6 +22,11 @@ pub fn tolerance<T: ComplexField>(n: usize) -> T::Real {
 pub fn close<T: ComplexField>(a: MatRef<'_, T>, b: MatRef<'_, T>) {
 	let error = (a - b).norm_l2();
 	let bound = tolerance::<T>(a.nrows()) * b.norm_l2().fmax(one::<T::Real>());
+	// An infinity on one side makes both `error` and `bound` infinite, and `inf <= inf` holds.
+	assert!(
+		error.is_finite() && bound.is_finite(),
+		"non-finite residual {error:?} or bound {bound:?}"
+	);
 	assert!(error <= bound, "residual {error:?}, bound {bound:?}");
 }
 pub fn unitary<T: ComplexField>(q: MatRef<'_, T>) {