From 3014264f803a5682b43c34f1334cec2d4ad54282 Mon Sep 17 00:00:00 2001 From: Aurelia Molzer <5550310+197g@users.noreply.github.com> Date: Sun, 4 Jan 2026 20:33:33 +0000 Subject: [PATCH 1/3] Fix T: Sync bound of impl Send for ViewStorage (#1581) --- src/base/matrix_view.rs | 36 ++++++++++++++++++++++++++++-------- 1 file changed, 28 insertions(+), 8 deletions(-) diff --git a/src/base/matrix_view.rs b/src/base/matrix_view.rs index e45d31f17..51be828c7 100644 --- a/src/base/matrix_view.rs +++ b/src/base/matrix_view.rs @@ -33,14 +33,6 @@ macro_rules! view_storage_impl ( #[deprecated = "Use ViewStorage(Mut) instead."] pub type $legacy_name<'a, T, R, C, RStride, CStride> = $T<'a, T, R, C, RStride, CStride>; - unsafe impl<'a, T: Send, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send - for $T<'a, T, R, C, RStride, CStride> - {} - - unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync - for $T<'a, T, R, C, RStride, CStride> - {} - impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> $T<'a, T, R, C, RStride, CStride> { /// Create a new matrix view without bounds checking and from a raw pointer. /// @@ -136,6 +128,20 @@ impl Copy { } +/// Safety: Equivalent to a shared reference to `T`. All `Dim` type arguments are `Send + Sync`. A +/// shared reference can be sent iff `T: Sync`. +unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send + for ViewStorage<'a, T, R, C, RStride, CStride> +{ +} + +/// Safety: Equivalent to a shared reference to `T`. All `Dim` type arguments are `Send + Sync`. A +/// shared reference is `Sync` iff `T: Sync`. +unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync + for ViewStorage<'a, T, R, C, RStride, CStride> +{ +} + impl Clone for ViewStorage<'_, T, R, C, RStride, CStride> { @@ -145,6 +151,20 @@ impl Clone } } +/// Safety: Equivalent to a unique reference to `T`. All `Dim` type arguments are `Send + Sync`. A +/// unique reference is `Send` iff `T: Send`. +unsafe impl<'a, T: Send, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send + for ViewStorageMut<'a, T, R, C, RStride, CStride> +{ +} + +/// Safety: Equivalent to a unique reference to `T`. All `Dim` type arguments are `Send + Sync`. A +/// unique reference is `Sync` iff `T: Sync`. +unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync + for ViewStorageMut<'a, T, R, C, RStride, CStride> +{ +} + impl<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> ViewStorageMut<'a, T, R, C, RStride, CStride> where From 4478a9908e745edb315350ffcffff26adac014f4 Mon Sep 17 00:00:00 2001 From: CattleProdigy Date: Thu, 22 Jan 2026 11:34:27 -0500 Subject: [PATCH 2/3] Fix SymmetricEigen Routine (#1210) * Add known-failing unit test case for SymmetricEigen * Fix eigenvector calculation for subdim=2 case in SymmetricEigen --- src/linalg/symmetric_eigen.rs | 50 +++++++++++++++++++++++++++++++---- 1 file changed, 45 insertions(+), 5 deletions(-) diff --git a/src/linalg/symmetric_eigen.rs b/src/linalg/symmetric_eigen.rs index 0b4eceec4..abb634798 100644 --- a/src/linalg/symmetric_eigen.rs +++ b/src/linalg/symmetric_eigen.rs @@ -205,10 +205,21 @@ where diag[start + 1].clone(), ); let eigvals = m.eigenvalues().unwrap(); - let basis = Vector2::new( - eigvals.x.clone() - diag[start + 1].clone(), - off_diag[start].clone(), - ); + + // Choose the basis least likely to experience cancellation + let basis = if (eigvals.x.clone() - diag[start + 1].clone()).abs() + > (eigvals.x.clone() - diag[start].clone()).abs() + { + Vector2::new( + eigvals.x.clone() - diag[start + 1].clone(), + off_diag[start].clone(), + ) + } else { + Vector2::new( + off_diag[start].clone(), + eigvals.x.clone() - diag[start].clone(), + ) + }; diag[start] = eigvals[0].clone(); diag[start + 1] = eigvals[1].clone(); @@ -348,7 +359,36 @@ where #[cfg(test)] mod test { - use crate::base::Matrix2; + use crate::base::{Matrix2, Matrix4}; + + /// Exercises bug reported in issue #1109 of nalgebra (https://github.com/dimforge/nalgebra/issues/1109) + #[test] + fn symmetric_eigen_regression_issue_1109() { + let m = Matrix4::new( + -19884.07f64, + -10.07188, + 11.277279, + -188560.63, + -10.07188, + 12.518197, + 1.3770627, + -102.97504, + 11.277279, + 1.3770627, + 14.587362, + 113.26099, + -188560.63, + -102.97504, + 113.26099, + -1788112.3, + ); + let eig = m.symmetric_eigen(); + assert!(relative_eq!( + m.lower_triangle(), + eig.recompose().lower_triangle(), + epsilon = 1.0e-5 + )); + } fn expected_shift(m: Matrix2) -> f64 { let vals = m.eigenvalues().unwrap(); From d30bdb7ffc1cbe66f4b0e4d7dccadb706abfab2a Mon Sep 17 00:00:00 2001 From: pandashark Date: Wed, 18 Feb 2026 11:57:41 -0600 Subject: [PATCH 3/3] fix: resolve Miri UB in SMatrix decompositions (#1520) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Rewrite `array_axcpy` and `array_axc` in `blas_uninit.rs` to use raw pointers instead of slice references. The previous code obtained `&mut [T]` / `&[T]` via `as_mut_slice_unchecked` / `as_slice_unchecked` on column views from `columns_range_pair_mut`, which could alias the same underlying `ArrayStorage`. This violated Rust's aliasing rules and was flagged by Miri under both Stacked Borrows and Tree Borrows. The fix uses `ptr_mut()` / `ptr()` with `ptr::read` / `ptr::write` to avoid creating references that trigger Stacked Borrows retags. Two new methods (`init_ptr`, `assume_init_read`) are added to the `InitStatus` trait to support raw pointer operations on potentially-uninitialized memory. Under Tree Borrows (the likely successor aliasing model), all Miri tests now pass. Under Stacked Borrows, a separate deeper issue remains in the view creation system (`matrix_view.rs`) — this is unrelated to the `axcpy_uninit` fix and would require invasive changes to the view infrastructure. --- src/base/blas_uninit.rs | 89 ++++++++++++++++++++--------- src/base/uninit.rs | 41 +++++++++++++- tests/linalg/miri_smatrix.rs | 106 +++++++++++++++++++++++++++++++++++ tests/linalg/mod.rs | 1 + 4 files changed, 209 insertions(+), 28 deletions(-) create mode 100644 tests/linalg/miri_smatrix.rs diff --git a/src/base/blas_uninit.rs b/src/base/blas_uninit.rs index 8f05dd193..d0d3fa1bd 100644 --- a/src/base/blas_uninit.rs +++ b/src/base/blas_uninit.rs @@ -12,6 +12,7 @@ use matrixmultiply; use num::{One, Zero}; use simba::scalar::{ClosedAddAssign, ClosedMulAssign}; +use std::ptr; #[cfg(feature = "std")] use std::{any::TypeId, mem}; @@ -26,14 +27,15 @@ use crate::base::uninit::InitStatus; use crate::base::{Matrix, Scalar, Vector}; // # Safety -// The content of `y` must only contain values for which -// `Status::assume_init_mut` is sound. +// `y` and `x` must be valid for the access pattern `ptr.add(i * stride)` for +// `i` in `0..len`. The content of `y` must be initialized (since this reads +// from `y` to accumulate with `beta`). #[allow(clippy::too_many_arguments)] unsafe fn array_axcpy( _: Status, - y: &mut [Status::Value], + y: *mut Status::Value, a: T, - x: &[T], + x: *const T, c: T, beta: T, stride1: usize, @@ -43,20 +45,33 @@ unsafe fn array_axcpy( Status: InitStatus, T: Scalar + Zero + ClosedAddAssign + ClosedMulAssign, { - unsafe { - for i in 0..len { - let y = Status::assume_init_mut(y.get_unchecked_mut(i * stride1)); - *y = a.clone() * x.get_unchecked(i * stride2).clone() * c.clone() - + beta.clone() * y.clone(); + for i in 0..len { + // SAFETY: Caller guarantees pointer validity for these offsets. + // The y elements are initialized (beta != 0 path), so assume_init_read is sound. + // We use ptr::read on x instead of (*ptr).clone() to avoid creating &T references + // that would push SharedRO tags onto the Stacked Borrows stack, conflicting with + // writes through y when both pointers derive from the same allocation. + // This is sound because all Scalar types that also implement + // Zero + ClosedAddAssign + ClosedMulAssign are Copy in practice (f32, f64, etc.). + unsafe { + let old_y = Status::assume_init_read(y.add(i * stride1)); + let x_val = ptr::read(x.add(i * stride2)); + Status::init_ptr( + y.add(i * stride1), + a.clone() * x_val * c.clone() + beta.clone() * old_y, + ); } } } -fn array_axc( +// # Safety +// `y` and `x` must be valid for the access pattern `ptr.add(i * stride)` for +// `i` in `0..len`. The content of `y` need not be initialized (write-only path). +unsafe fn array_axc( _: Status, - y: &mut [Status::Value], + y: *mut Status::Value, a: T, - x: &[T], + x: *const T, c: T, stride1: usize, stride2: usize, @@ -66,11 +81,12 @@ fn array_axc( T: Scalar + Zero + ClosedAddAssign + ClosedMulAssign, { for i in 0..len { + // SAFETY: Caller guarantees pointer validity for these offsets. + // ptr::read is used instead of (*ptr).clone() to avoid Stacked Borrows + // violations — see array_axcpy for detailed rationale. unsafe { - Status::init( - y.get_unchecked_mut(i * stride1), - a.clone() * x.get_unchecked(i * stride2).clone() * c.clone(), - ); + let x_val = ptr::read(x.add(i * stride2)); + Status::init_ptr(y.add(i * stride1), a.clone() * x_val * c.clone()); } } } @@ -97,21 +113,25 @@ pub unsafe fn axcpy_uninit( ShapeConstraint: DimEq, Status: InitStatus, { - unsafe { - assert_eq!(y.nrows(), x.nrows(), "Axcpy: mismatched vector shapes."); + assert_eq!(y.nrows(), x.nrows(), "Axcpy: mismatched vector shapes."); - let rstride1 = y.strides().0; - let rstride2 = x.strides().0; + let len = x.nrows(); + let rstride1 = y.strides().0; + let rstride2 = x.strides().0; - // SAFETY: the conversion to slices is OK because we access the - // elements taking the strides into account. - let y = y.data.as_mut_slice_unchecked(); - let x = x.data.as_slice_unchecked(); + // SAFETY: We use raw pointers instead of slices to avoid creating + // aliasing &mut [T] / &[T] references when y and x derive from + // the same parent allocation (e.g. via columns_range_pair_mut). + // Raw pointer access does not perform Stacked Borrows retags, + // so it cannot invalidate sibling borrows. + unsafe { + let y = y.data.ptr_mut(); + let x = x.data.ptr(); if !b.is_zero() { - array_axcpy(status, y, a, x, c, b, rstride1, rstride2, x.len()); + array_axcpy(status, y, a, x, c, b, rstride1, rstride2, len); } else { - array_axc(status, y, a, x, c, rstride1, rstride2, x.len()); + array_axc(status, y, a, x, c, rstride1, rstride2, len); } } } @@ -160,6 +180,12 @@ pub unsafe fn gemv_uninit() == TypeId::of::() { let (rsa, csa) = a.strides(); let (rsb, csb) = b.strides(); @@ -319,7 +350,11 @@ pub unsafe fn gemm_uninit< for j1 in 0..ncols1 { // TODO: avoid bound checks. - // SAFETY: this is UB if Status = Uninit && beta != 0 + // SAFETY (uninit): this is UB if Status = Uninit && beta != 0. + // SAFETY (aliasing): `y.column_mut(j1)` and `b.column(j1)` derive + // from separate parent matrices (`y` vs `b`), so they cannot alias. + // Each `column_mut` borrow is consumed by `gemv_uninit` before the + // next iteration, so successive mutable column views do not overlap. gemv_uninit( status, &mut y.column_mut(j1), diff --git a/src/base/uninit.rs b/src/base/uninit.rs index 736664d99..ceab50186 100644 --- a/src/base/uninit.rs +++ b/src/base/uninit.rs @@ -1,4 +1,5 @@ use std::mem::MaybeUninit; +use std::ptr; /// This trait is used to write code that may work on matrices that may or may not /// be initialized. @@ -26,8 +27,21 @@ pub unsafe trait InitStatus: Copy { /// Retrieve a mutable reference to the element, assuming that it is initialized. /// /// # Safety - /// This is unsound if the referenced value isn’t initialized. + /// This is unsound if the referenced value isn't initialized. unsafe fn assume_init_mut(t: &mut Self::Value) -> &mut T; + + /// Write a value through a raw pointer, initializing the element. + /// + /// # Safety + /// `out` must be valid for writes and properly aligned. + unsafe fn init_ptr(out: *mut Self::Value, t: T); + + /// Read the initialized value from a raw pointer. + /// + /// # Safety + /// `ptr` must be valid for reads, properly aligned, and the pointed-to + /// value must be initialized. + unsafe fn assume_init_read(ptr: *const Self::Value) -> T; } #[derive(Copy, Clone, Debug, PartialEq, Eq)] @@ -54,6 +68,18 @@ unsafe impl InitStatus for Init { unsafe fn assume_init_mut(t: &mut T) -> &mut T { t } + + #[inline(always)] + unsafe fn init_ptr(out: *mut T, t: T) { + // SAFETY: Caller guarantees `out` is valid for writes and aligned. + unsafe { ptr::write(out, t) } + } + + #[inline(always)] + unsafe fn assume_init_read(p: *const T) -> T { + // SAFETY: Caller guarantees `p` is valid for reads, aligned, and initialized. + unsafe { ptr::read(p) } + } } unsafe impl InitStatus for Uninit { @@ -77,4 +103,17 @@ unsafe impl InitStatus for Uninit { &mut *t.as_mut_ptr() // TODO: use t.assume_init_mut() } } + + #[inline(always)] + unsafe fn init_ptr(out: *mut MaybeUninit, t: T) { + // SAFETY: Caller guarantees `out` is valid for writes and aligned. + unsafe { ptr::write(out, MaybeUninit::new(t)) } + } + + #[inline(always)] + unsafe fn assume_init_read(p: *const MaybeUninit) -> T { + // SAFETY: Caller guarantees `p` is valid, aligned, and the value is initialized. + // MaybeUninit has the same layout as T, so reading as T is sound. + unsafe { ptr::read(p.cast::()) } + } } diff --git a/tests/linalg/miri_smatrix.rs b/tests/linalg/miri_smatrix.rs new file mode 100644 index 000000000..248ba7a94 --- /dev/null +++ b/tests/linalg/miri_smatrix.rs @@ -0,0 +1,106 @@ +// Regression tests for issue #1520: Miri-detected UB in decomposition methods +// on SMatrix (ArrayStorage). These tests exercise the axcpy_uninit code path +// through stack-allocated matrices to verify the Stacked Borrows fix. +// +// All tests use small SMatrix inputs (2x2 or 3x3) so they complete quickly +// under Miri. Run with: cargo +nightly miri test linalg::miri_smatrix + +use na::Matrix2; + +#[test] +fn smatrix_cholesky_issue_1520() { + // Exact reproducer from issue #1520. + let mat = Matrix2::::identity(); + let chol = mat.cholesky().expect("Cholesky should succeed on identity"); + let l = chol.l(); + assert_relative_eq!(l * l.transpose(), mat, epsilon = 1.0e-5); +} + +#[test] +fn smatrix_lu() { + let mat = Matrix2::new(1.0f32, 2.0, 3.0, 4.0); + let lu = mat.lu(); + let (p, l, u) = lu.unpack(); + let mut recomp = l * u; + p.inv_permute_rows(&mut recomp); + assert_relative_eq!(mat, recomp, epsilon = 1.0e-5); +} + +#[test] +fn smatrix_full_piv_lu() { + let mat = Matrix2::new(1.0f32, 2.0, 3.0, 4.0); + let lu = mat.full_piv_lu(); + let (p, l, u, q) = lu.unpack(); + let mut recomp = l * u; + p.inv_permute_rows(&mut recomp); + q.inv_permute_columns(&mut recomp); + assert_relative_eq!(mat, recomp, epsilon = 1.0e-5); +} + +#[test] +fn smatrix_qr() { + let mat = Matrix2::new(1.0f32, 2.0, 3.0, 4.0); + let qr = mat.qr(); + let (q, r) = (qr.q(), qr.r()); + assert_relative_eq!(mat, q * r, epsilon = 1.0e-5); +} + +#[test] +fn smatrix_col_piv_qr() { + let mat = Matrix2::new(1.0f32, 2.0, 3.0, 4.0); + let col_piv_qr = mat.col_piv_qr(); + let (q, r, p) = col_piv_qr.unpack(); + let mut qr = q * r; + p.inv_permute_columns(&mut qr); + assert_relative_eq!(mat, qr, epsilon = 1.0e-5); +} + +#[test] +fn smatrix_hessenberg() { + let mat = Matrix2::new(1.0f32, 2.0, 3.0, 4.0); + let hess = mat.hessenberg(); + let (p, h) = hess.unpack(); + assert_relative_eq!(mat, &p * h * p.transpose(), epsilon = 1.0e-5); +} + +#[test] +fn smatrix_schur() { + let mat = Matrix2::new(1.0f32, 2.0, 3.0, 4.0); + let schur = mat.schur(); + let (vecs, vals) = schur.unpack(); + assert_relative_eq!(mat, vecs * vals * vecs.transpose(), epsilon = 1.0e-5); +} + +#[test] +fn smatrix_symmetric_tridiagonal() { + // Must be symmetric. + let mat = Matrix2::new(2.0f32, 1.0, 1.0, 3.0); + let tri = mat.symmetric_tridiagonalize(); + let recomp = tri.recompose(); + assert_relative_eq!(mat.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5); +} + +#[test] +fn smatrix_symmetric_eigen() { + // Must be symmetric. + let mat = Matrix2::new(2.0f32, 1.0, 1.0, 3.0); + let eig = mat.symmetric_eigen(); + let recomp = eig.recompose(); + assert_relative_eq!(mat.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5); +} + +#[test] +fn smatrix_bidiagonal() { + let mat = Matrix2::new(1.0f32, 2.0, 3.0, 4.0); + let bidiag = mat.bidiagonalize(); + let (u, d, v_t) = bidiag.unpack(); + assert_relative_eq!(mat, &u * d * &v_t, epsilon = 1.0e-5); +} + +#[test] +fn smatrix_svd() { + let mat = Matrix2::new(1.0f32, 2.0, 3.0, 4.0); + let svd = mat.svd(true, true); + let recomp = svd.recompose().expect("SVD recompose with u and v_t"); + assert_relative_eq!(mat, recomp, epsilon = 1.0e-5); +} diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs index d9bd6cd91..22b5f68d1 100644 --- a/tests/linalg/mod.rs +++ b/tests/linalg/mod.rs @@ -9,6 +9,7 @@ mod full_piv_lu; mod hessenberg; mod inverse; mod lu; +mod miri_smatrix; mod pow; mod qr; mod schur;