diff --git a/CHANGELOG.md b/CHANGELOG.md index 04ea1c34..8eae0834 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,13 @@ documented here. This project adheres to [Semantic Versioning](https://semver.org/). +## [0.29.0] - WIP +### Modified +- The closure given to `apply`, `zip_apply`, `zip_zip_apply` must now modify the + first argument inplace, instead of returning a new value. This makes these + methods more versatile, and avoid useless clones when using non-Copy scalar + types. + ## [0.28.0] ### Added - Implement `Hash` for `Transform`. diff --git a/Cargo.toml b/Cargo.toml index d10db84a..9c433b2a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -31,7 +31,6 @@ io = [ "pest", "pest_derive" ] compare = [ "matrixcompare-core" ] libm = [ "simba/libm" ] libm-force = [ "simba/libm_force" ] -no_unsound_assume_init = [ ] macros = [ "nalgebra-macros" ] # Conversion diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 1bca79a5..f6628bfe 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -9,7 +9,6 @@ use simba::scalar::RealField; use crate::ComplexHelper; use na::allocator::Allocator; use na::dimension::{Const, Dim}; -use na::storage::Storage; use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; use lapack; @@ -73,14 +72,15 @@ where let ljob = if left_eigenvectors { b'V' } else { b'T' }; let rjob = if eigenvectors { b'V' } else { b'T' }; - let (nrows, ncols) = m.data.shape(); + let (nrows, ncols) = m.shape_generic(); let n = nrows.value(); let lda = n as i32; - let mut wr = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; + // TODO: avoid the initialization? + let mut wr = Matrix::zeros_generic(nrows, Const::<1>); // TODO: Tap into the workspace. - let mut wi = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; + let mut wi = Matrix::zeros_generic(nrows, Const::<1>); let mut info = 0; let mut placeholder1 = [T::zero()]; @@ -103,14 +103,13 @@ where lapack_check!(info); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; match (left_eigenvectors, eigenvectors) { (true, true) => { - let mut vl = - unsafe { Matrix::new_uninitialized_generic(nrows, ncols).assume_init() }; - let mut vr = - unsafe { Matrix::new_uninitialized_generic(nrows, ncols).assume_init() }; + // TODO: avoid the initializations? + let mut vl = Matrix::zeros_generic(nrows, ncols); + let mut vr = Matrix::zeros_generic(nrows, ncols); T::xgeev( ljob, @@ -139,8 +138,8 @@ where } } (true, false) => { - let mut vl = - unsafe { Matrix::new_uninitialized_generic(nrows, ncols).assume_init() }; + // TODO: avoid the initialization? + let mut vl = Matrix::zeros_generic(nrows, ncols); T::xgeev( ljob, @@ -169,8 +168,8 @@ where } } (false, true) => { - let mut vr = - unsafe { Matrix::new_uninitialized_generic(nrows, ncols).assume_init() }; + // TODO: avoid the initialization? + let mut vr = Matrix::zeros_generic(nrows, ncols); T::xgeev( ljob, @@ -242,13 +241,14 @@ where "Unable to compute the eigenvalue decomposition of a non-square matrix." ); - let nrows = m.data.shape().0; + let nrows = m.shape_generic().0; let n = nrows.value(); let lda = n as i32; - let mut wr = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; - let mut wi = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; + // TODO: avoid the initialization? + let mut wr = Matrix::zeros_generic(nrows, Const::<1>); + let mut wi = Matrix::zeros_generic(nrows, Const::<1>); let mut info = 0; let mut placeholder1 = [T::zero()]; @@ -271,7 +271,7 @@ where lapack_panic!(info); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; T::xgeev( b'T', @@ -291,7 +291,7 @@ where ); lapack_panic!(info); - let mut res = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; + let mut res = Matrix::zeros_generic(nrows, Const::<1>); for i in 0..res.len() { res[i] = Complex::new(wr[i], wi[i]); diff --git a/nalgebra-lapack/src/hessenberg.rs b/nalgebra-lapack/src/hessenberg.rs index 0d911cf8..e05349d9 100644 --- a/nalgebra-lapack/src/hessenberg.rs +++ b/nalgebra-lapack/src/hessenberg.rs @@ -4,7 +4,6 @@ use num_complex::Complex; use crate::ComplexHelper; use na::allocator::Allocator; use na::dimension::{Const, DimDiff, DimSub, U1}; -use na::storage::Storage; use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; use lapack; @@ -48,7 +47,7 @@ where { /// Computes the hessenberg decomposition of the matrix `m`. pub fn new(mut m: OMatrix) -> Self { - let nrows = m.data.shape().0; + let nrows = m.shape_generic().0; let n = nrows.value() as i32; assert!( @@ -60,14 +59,12 @@ where "Unable to compute the hessenberg decomposition of an empty matrix." ); - let mut tau = unsafe { - Matrix::new_uninitialized_generic(nrows.sub(Const::<1>), Const::<1>).assume_init() - }; + let mut tau = Matrix::zeros_generic(nrows.sub(Const::<1>), Const::<1>); let mut info = 0; let lwork = T::xgehrd_work_size(n, 1, n, m.as_mut_slice(), n, tau.as_mut_slice(), &mut info); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; lapack_panic!(info); diff --git a/nalgebra-lapack/src/lib.rs b/nalgebra-lapack/src/lib.rs index 9a027772..84fa03fa 100644 --- a/nalgebra-lapack/src/lib.rs +++ b/nalgebra-lapack/src/lib.rs @@ -139,10 +139,3 @@ impl ComplexHelper for Complex { self.re } } - -unsafe fn uninitialized_vec(n: usize) -> Vec { - let mut res = Vec::new(); - res.reserve_exact(n); - res.set_len(n); - res -} diff --git a/nalgebra-lapack/src/lu.rs b/nalgebra-lapack/src/lu.rs index 2130fc7e..7540c75e 100644 --- a/nalgebra-lapack/src/lu.rs +++ b/nalgebra-lapack/src/lu.rs @@ -61,7 +61,7 @@ where { /// Computes the LU decomposition with partial (row) pivoting of `matrix`. pub fn new(mut m: OMatrix) -> Self { - let (nrows, ncols) = m.data.shape(); + let (nrows, ncols) = m.shape_generic(); let min_nrows_ncols = nrows.min(ncols); let nrows = nrows.value() as i32; let ncols = ncols.value() as i32; @@ -87,7 +87,7 @@ where #[inline] #[must_use] pub fn l(&self) -> OMatrix> { - let (nrows, ncols) = self.lu.data.shape(); + let (nrows, ncols) = self.lu.shape_generic(); let mut res = self.lu.columns_generic(0, nrows.min(ncols)).into_owned(); res.fill_upper_triangle(Zero::zero(), 1); @@ -100,7 +100,7 @@ where #[inline] #[must_use] pub fn u(&self) -> OMatrix, C> { - let (nrows, ncols) = self.lu.data.shape(); + let (nrows, ncols) = self.lu.shape_generic(); let mut res = self.lu.rows_generic(0, nrows.min(ncols)).into_owned(); res.fill_lower_triangle(Zero::zero(), 1); @@ -115,7 +115,7 @@ where #[inline] #[must_use] pub fn p(&self) -> OMatrix { - let (dim, _) = self.lu.data.shape(); + let (dim, _) = self.lu.shape_generic(); let mut id = Matrix::identity_generic(dim, dim); self.permute(&mut id); @@ -290,7 +290,7 @@ where ); lapack_check!(info); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; T::xgetri( dim, diff --git a/nalgebra-lapack/src/qr.rs b/nalgebra-lapack/src/qr.rs index 7f00d058..895e34f3 100644 --- a/nalgebra-lapack/src/qr.rs +++ b/nalgebra-lapack/src/qr.rs @@ -7,7 +7,6 @@ use num_complex::Complex; use crate::ComplexHelper; use na::allocator::Allocator; use na::dimension::{Const, Dim, DimMin, DimMinimum}; -use na::storage::Storage; use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; use lapack; @@ -54,12 +53,10 @@ where { /// Computes the QR decomposition of the matrix `m`. pub fn new(mut m: OMatrix) -> Self { - let (nrows, ncols) = m.data.shape(); + let (nrows, ncols) = m.shape_generic(); let mut info = 0; - let mut tau = unsafe { - Matrix::new_uninitialized_generic(nrows.min(ncols), Const::<1>).assume_init() - }; + let mut tau = Matrix::zeros_generic(nrows.min(ncols), Const::<1>); if nrows.value() == 0 || ncols.value() == 0 { return Self { qr: m, tau }; @@ -74,7 +71,7 @@ where &mut info, ); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; T::xgeqrf( nrows.value() as i32, @@ -94,7 +91,7 @@ where #[inline] #[must_use] pub fn r(&self) -> OMatrix, C> { - let (nrows, ncols) = self.qr.data.shape(); + let (nrows, ncols) = self.qr.shape_generic(); self.qr.rows_generic(0, nrows.min(ncols)).upper_triangle() } } @@ -120,7 +117,7 @@ where #[inline] #[must_use] pub fn q(&self) -> OMatrix> { - let (nrows, ncols) = self.qr.data.shape(); + let (nrows, ncols) = self.qr.shape_generic(); let min_nrows_ncols = nrows.min(ncols); if min_nrows_ncols.value() == 0 { diff --git a/nalgebra-lapack/src/schur.rs b/nalgebra-lapack/src/schur.rs index 644f8a5c..13dfc05e 100644 --- a/nalgebra-lapack/src/schur.rs +++ b/nalgebra-lapack/src/schur.rs @@ -9,7 +9,6 @@ use simba::scalar::RealField; use crate::ComplexHelper; use na::allocator::Allocator; use na::dimension::{Const, Dim}; -use na::storage::Storage; use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; use lapack; @@ -71,16 +70,16 @@ where "Unable to compute the eigenvalue decomposition of a non-square matrix." ); - let (nrows, ncols) = m.data.shape(); + let (nrows, ncols) = m.shape_generic(); let n = nrows.value(); let lda = n as i32; let mut info = 0; - let mut wr = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; - let mut wi = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; - let mut q = unsafe { Matrix::new_uninitialized_generic(nrows, ncols).assume_init() }; + let mut wr = Matrix::zeros_generic(nrows, Const::<1>); + let mut wi = Matrix::zeros_generic(nrows, Const::<1>); + let mut q = Matrix::zeros_generic(nrows, ncols); // Placeholders: let mut bwork = [0i32]; let mut unused = 0; @@ -101,7 +100,7 @@ where ); lapack_check!(info); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; T::xgees( b'V', @@ -153,9 +152,7 @@ where where DefaultAllocator: Allocator, D>, { - let mut out = unsafe { - OVector::new_uninitialized_generic(self.t.data.shape().0, Const::<1>).assume_init() - }; + let mut out = Matrix::zeros_generic(self.t.shape_generic().0, Const::<1>); for i in 0..out.len() { out[i] = Complex::new(self.re[i], self.im[i]) diff --git a/nalgebra-lapack/src/svd.rs b/nalgebra-lapack/src/svd.rs index 3357e621..972ffa1b 100644 --- a/nalgebra-lapack/src/svd.rs +++ b/nalgebra-lapack/src/svd.rs @@ -6,7 +6,6 @@ use std::cmp; use na::allocator::Allocator; use na::dimension::{Const, Dim, DimMin, DimMinimum, U1}; -use na::storage::Storage; use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; use lapack; @@ -89,7 +88,7 @@ macro_rules! svd_impl( Allocator<$t, DimMinimum> { fn compute(mut m: OMatrix<$t, R, C>) -> Option> { - let (nrows, ncols) = m.data.shape(); + let (nrows, ncols) = m.shape_generic(); if nrows.value() == 0 || ncols.value() == 0 { return None; @@ -99,9 +98,9 @@ macro_rules! svd_impl( let lda = nrows.value() as i32; - let mut u = unsafe { Matrix::new_uninitialized_generic(nrows, nrows).assume_init() }; - let mut s = unsafe { Matrix::new_uninitialized_generic(nrows.min(ncols), Const::<1>).assume_init() }; - let mut vt = unsafe { Matrix::new_uninitialized_generic(ncols, ncols).assume_init() }; + let mut u = Matrix::zeros_generic(nrows, nrows); + let mut s = Matrix::zeros_generic(nrows.min(ncols), Const::<1>); + let mut vt = Matrix::zeros_generic(ncols, ncols); let ldu = nrows.value(); let ldvt = ncols.value(); @@ -109,7 +108,7 @@ macro_rules! svd_impl( let mut work = [ 0.0 ]; let mut lwork = -1 as i32; let mut info = 0; - let mut iwork = unsafe { crate::uninitialized_vec(8 * cmp::min(nrows.value(), ncols.value())) }; + let mut iwork = vec![0; 8 * cmp::min(nrows.value(), ncols.value())]; unsafe { $lapack_func(job, nrows.value() as i32, ncols.value() as i32, m.as_mut_slice(), @@ -119,7 +118,7 @@ macro_rules! svd_impl( lapack_check!(info); lwork = work[0] as i32; - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![0.0; lwork as usize]; unsafe { $lapack_func(job, nrows.value() as i32, ncols.value() as i32, m.as_mut_slice(), @@ -151,8 +150,8 @@ macro_rules! svd_impl( /// been manually changed by the user. #[inline] pub fn recompose(self) -> OMatrix<$t, R, C> { - let nrows = self.u.data.shape().0; - let ncols = self.vt.data.shape().1; + let nrows = self.u.shape_generic().0; + let ncols = self.vt.shape_generic().1; let min_nrows_ncols = nrows.min(ncols); let mut res: OMatrix<_, R, C> = Matrix::zeros_generic(nrows, ncols); @@ -177,8 +176,8 @@ macro_rules! svd_impl( #[inline] #[must_use] pub fn pseudo_inverse(&self, epsilon: $t) -> OMatrix<$t, C, R> { - let nrows = self.u.data.shape().0; - let ncols = self.vt.data.shape().1; + let nrows = self.u.shape_generic().0; + let ncols = self.vt.shape_generic().1; let min_nrows_ncols = nrows.min(ncols); let mut res: OMatrix<_, C, R> = Matrix::zeros_generic(ncols, nrows); @@ -241,7 +240,7 @@ macro_rules! svd_complex_impl( Allocator, R, R> + Allocator, C, C> + Allocator<$t, DimMinimum> { - let (nrows, ncols) = m.data.shape(); + let (nrows, ncols) = m.shape_generic(); if nrows.value() == 0 || ncols.value() == 0 { return None; @@ -254,9 +253,9 @@ macro_rules! svd_complex_impl( let min_nrows_ncols = nrows.min(ncols); - let mut u = unsafe { Matrix::new_uninitialized_generic(nrows, nrows) }; - let mut s = unsafe { Matrix::new_uninitialized_generic(min_nrows_ncols, U1) }; - let mut vt = unsafe { Matrix::new_uninitialized_generic(ncols, ncols) }; + let mut u = Matrix::zeros_generic(nrows, nrows); + let mut s = Matrix::zeros_generic(min_nrows_ncols, U1); + let mut vt = Matrix::zeros_generic(ncols, ncols); let ldu = nrows.value(); let ldvt = ncols.value(); diff --git a/nalgebra-lapack/src/symmetric_eigen.rs b/nalgebra-lapack/src/symmetric_eigen.rs index d276437e..8cbe63f8 100644 --- a/nalgebra-lapack/src/symmetric_eigen.rs +++ b/nalgebra-lapack/src/symmetric_eigen.rs @@ -9,7 +9,6 @@ use simba::scalar::RealField; use crate::ComplexHelper; use na::allocator::Allocator; use na::dimension::{Const, Dim}; -use na::storage::Storage; use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; use lapack; @@ -89,19 +88,18 @@ where let jobz = if eigenvectors { b'V' } else { b'T' }; - let nrows = m.data.shape().0; + let nrows = m.shape_generic().0; let n = nrows.value(); let lda = n as i32; - let mut values = - unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; + let mut values = Matrix::zeros_generic(nrows, Const::<1>); let mut info = 0; let lwork = T::xsyev_work_size(jobz, b'L', n as i32, m.as_mut_slice(), lda, &mut info); lapack_check!(info); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; T::xsyev( jobz, diff --git a/nalgebra-sparse/src/convert/impl_std_ops.rs b/nalgebra-sparse/src/convert/impl_std_ops.rs index ba4c015b..683227e2 100644 --- a/nalgebra-sparse/src/convert/impl_std_ops.rs +++ b/nalgebra-sparse/src/convert/impl_std_ops.rs @@ -2,7 +2,7 @@ use crate::convert::serial::*; use crate::coo::CooMatrix; use crate::csc::CscMatrix; use crate::csr::CsrMatrix; -use nalgebra::storage::Storage; +use nalgebra::storage::RawStorage; use nalgebra::{ClosedAdd, DMatrix, Dim, Matrix, Scalar}; use num_traits::Zero; @@ -11,7 +11,7 @@ where T: Scalar + Zero, R: Dim, C: Dim, - S: Storage, + S: RawStorage, { fn from(matrix: &'a Matrix) -> Self { convert_dense_coo(matrix) @@ -50,7 +50,7 @@ where T: Scalar + Zero, R: Dim, C: Dim, - S: Storage, + S: RawStorage, { fn from(matrix: &'a Matrix) -> Self { convert_dense_csr(matrix) @@ -89,7 +89,7 @@ where T: Scalar + Zero, R: Dim, C: Dim, - S: Storage, + S: RawStorage, { fn from(matrix: &'a Matrix) -> Self { convert_dense_csc(matrix) diff --git a/nalgebra-sparse/src/convert/serial.rs b/nalgebra-sparse/src/convert/serial.rs index 7e0da7bc..f84a6583 100644 --- a/nalgebra-sparse/src/convert/serial.rs +++ b/nalgebra-sparse/src/convert/serial.rs @@ -7,7 +7,7 @@ use std::ops::Add; use num_traits::Zero; -use nalgebra::storage::Storage; +use nalgebra::storage::RawStorage; use nalgebra::{ClosedAdd, DMatrix, Dim, Matrix, Scalar}; use crate::coo::CooMatrix; @@ -21,7 +21,7 @@ where T: Scalar + Zero, R: Dim, C: Dim, - S: Storage, + S: RawStorage, { let mut coo = CooMatrix::new(dense.nrows(), dense.ncols()); @@ -96,7 +96,7 @@ where T: Scalar + Zero, R: Dim, C: Dim, - S: Storage, + S: RawStorage, { let mut row_offsets = Vec::with_capacity(dense.nrows() + 1); let mut col_idx = Vec::new(); @@ -173,7 +173,7 @@ where T: Scalar + Zero, R: Dim, C: Dim, - S: Storage, + S: RawStorage, { let mut col_offsets = Vec::with_capacity(dense.ncols() + 1); let mut row_idx = Vec::new(); diff --git a/nalgebra-sparse/src/coo.rs b/nalgebra-sparse/src/coo.rs index 679dbdb2..34e5ceec 100644 --- a/nalgebra-sparse/src/coo.rs +++ b/nalgebra-sparse/src/coo.rs @@ -57,7 +57,7 @@ impl CooMatrix { /// Panics if any part of the dense matrix is out of bounds of the sparse matrix /// when inserted at `(r, c)`. #[inline] - pub fn push_matrix>( + pub fn push_matrix>( &mut self, r: usize, c: usize, diff --git a/nalgebra-sparse/src/ops/impl_std_ops.rs b/nalgebra-sparse/src/ops/impl_std_ops.rs index 590bd934..721023a5 100644 --- a/nalgebra-sparse/src/ops/impl_std_ops.rs +++ b/nalgebra-sparse/src/ops/impl_std_ops.rs @@ -7,7 +7,7 @@ use crate::ops::serial::{ }; use crate::ops::Op; use nalgebra::allocator::Allocator; -use nalgebra::base::storage::Storage; +use nalgebra::base::storage::RawStorage; use nalgebra::constraint::{DimEq, ShapeConstraint}; use nalgebra::{ ClosedAdd, ClosedDiv, ClosedMul, ClosedSub, DefaultAllocator, Dim, Dynamic, Matrix, OMatrix, @@ -272,7 +272,7 @@ macro_rules! impl_spmm_cs_dense { ($matrix_type_name:ident, $spmm_fn:ident) => { // Implement ref-ref impl_spmm_cs_dense!(&'a $matrix_type_name, &'a Matrix, $spmm_fn, |lhs, rhs| { - let (_, ncols) = rhs.data.shape(); + let (_, ncols) = rhs.shape_generic(); let nrows = Dynamic::new(lhs.nrows()); let mut result = OMatrix::::zeros_generic(nrows, ncols); $spmm_fn(T::zero(), &mut result, T::one(), Op::NoOp(lhs), Op::NoOp(rhs)); @@ -301,14 +301,14 @@ macro_rules! impl_spmm_cs_dense { T: Scalar + ClosedMul + ClosedAdd + ClosedSub + ClosedDiv + Neg + Zero + One, R: Dim, C: Dim, - S: Storage, + S: RawStorage, DefaultAllocator: Allocator, // TODO: Is it possible to simplify these bounds? ShapeConstraint: // Bounds so that we can turn OMatrix into a DMatrixSliceMut - DimEq>::Buffer as Storage>::RStride> + DimEq>::Buffer as RawStorage>::RStride> + DimEq - + DimEq>::Buffer as Storage>::CStride> + + DimEq>::Buffer as RawStorage>::CStride> // Bounds so that we can turn &Matrix into a DMatrixSlice + DimEq + DimEq diff --git a/nalgebra-sparse/src/ops/serial/mod.rs b/nalgebra-sparse/src/ops/serial/mod.rs index 87285525..d8f1a343 100644 --- a/nalgebra-sparse/src/ops/serial/mod.rs +++ b/nalgebra-sparse/src/ops/serial/mod.rs @@ -8,7 +8,6 @@ //! some operations which will be able to dynamically adapt the output pattern to fit the //! result, but these have yet to be implemented. -#[macro_use] macro_rules! assert_compatible_spmm_dims { ($c:expr, $a:expr, $b:expr) => {{ use crate::ops::Op::{NoOp, Transpose}; @@ -37,7 +36,6 @@ macro_rules! assert_compatible_spmm_dims { }}; } -#[macro_use] macro_rules! assert_compatible_spadd_dims { ($c:expr, $a:expr) => { use crate::ops::Op; diff --git a/src/base/alias.rs b/src/base/alias.rs index 6bc04813..68829d9a 100644 --- a/src/base/alias.rs +++ b/src/base/alias.rs @@ -5,6 +5,8 @@ use crate::base::storage::Owned; #[cfg(any(feature = "std", feature = "alloc"))] use crate::base::vec_storage::VecStorage; use crate::base::{ArrayStorage, Const, Matrix, Unit}; +use crate::storage::OwnedUninit; +use std::mem::MaybeUninit; /* * @@ -19,6 +21,9 @@ use crate::base::{ArrayStorage, Const, Matrix, Unit}; /// **Because this is an alias, not all its methods are listed here. See the [`Matrix`](crate::base::Matrix) type too.** pub type OMatrix = Matrix>; +/// An owned matrix with uninitialized data. +pub type UninitMatrix = Matrix, R, C, OwnedUninit>; + /// An owned matrix column-major matrix with `R` rows and `C` columns. /// /// **Because this is an alias, not all its methods are listed here. See the [`Matrix`](crate::base::Matrix) type too.** @@ -278,6 +283,9 @@ pub type OVector = Matrix>; /// A statically sized D-dimensional column vector. pub type SVector = Matrix, U1, ArrayStorage>; // Owned, U1>>; +/// An owned matrix with uninitialized data. +pub type UninitVector = Matrix, D, U1, OwnedUninit>; + /// An owned matrix column-major matrix with `R` rows and `C` columns. /// /// **Because this is an alias, not all its methods are listed here. See the [`Matrix`](crate::base::Matrix) type too.** diff --git a/src/base/allocator.rs b/src/base/allocator.rs index b0f6537b..29286420 100644 --- a/src/base/allocator.rs +++ b/src/base/allocator.rs @@ -1,12 +1,14 @@ //! Abstract definition of a matrix data storage allocator. use std::any::Any; -use std::mem; use crate::base::constraint::{SameNumberOfColumns, SameNumberOfRows, ShapeConstraint}; use crate::base::dimension::{Dim, U1}; -use crate::base::storage::ContiguousStorageMut; use crate::base::{DefaultAllocator, Scalar}; +use crate::storage::{IsContiguous, RawStorageMut}; +use crate::StorageMut; +use std::fmt::Debug; +use std::mem::MaybeUninit; /// A matrix allocator of a memory buffer that may contain `R::to_usize() * C::to_usize()` /// elements of type `T`. @@ -17,12 +19,21 @@ use crate::base::{DefaultAllocator, Scalar}; /// /// Every allocator must be both static and dynamic. Though not all implementations may share the /// same `Buffer` type. -pub trait Allocator: Any + Sized { +pub trait Allocator: Any + Sized { /// The type of buffer this allocator can instanciate. - type Buffer: ContiguousStorageMut + Clone; + type Buffer: StorageMut + IsContiguous + Clone + Debug; + /// The type of buffer with uninitialized components this allocator can instanciate. + type BufferUninit: RawStorageMut, R, C> + IsContiguous; /// Allocates a buffer with the given number of rows and columns without initializing its content. - unsafe fn allocate_uninitialized(nrows: R, ncols: C) -> mem::MaybeUninit; + fn allocate_uninit(nrows: R, ncols: C) -> Self::BufferUninit; + + /// Assumes a data buffer to be initialized. + /// + /// # Safety + /// The user must make sure that every single entry of the buffer has been initialized, + /// or Undefined Behavior will immediately occur. + unsafe fn assume_init(uninit: Self::BufferUninit) -> Self::Buffer; /// Allocates a buffer initialized with the content of the given iterator. fn allocate_from_iterator>( @@ -41,15 +52,15 @@ pub trait Reallocator: /// `buf`. Data stored by `buf` are linearly copied to the output: /// /// # Safety - /// * The copy is performed as if both were just arrays (without a matrix structure). - /// * If `buf` is larger than the output size, then extra elements of `buf` are truncated. - /// * If `buf` is smaller than the output size, then extra elements of the output are left - /// uninitialized. + /// The following invariants must be respected by the implementors of this method: + /// * The copy is performed as if both were just arrays (without taking into account the matrix structure). + /// * If the underlying buffer is being shrunk, the removed elements must **not** be dropped + /// by this method. Dropping them is the responsibility of the caller. unsafe fn reallocate_copy( nrows: RTo, ncols: CTo, buf: >::Buffer, - ) -> >::Buffer; + ) -> >::BufferUninit; } /// The number of rows of the result of a componentwise operation on two matrices. @@ -67,7 +78,6 @@ where R2: Dim, C1: Dim, C2: Dim, - T: Scalar, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { } @@ -78,7 +88,6 @@ where R2: Dim, C1: Dim, C2: Dim, - T: Scalar, DefaultAllocator: Allocator + Allocator, SameShapeC>, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { @@ -91,7 +100,6 @@ pub trait SameShapeVectorAllocator: where R1: Dim, R2: Dim, - T: Scalar, ShapeConstraint: SameNumberOfRows, { } @@ -100,7 +108,6 @@ impl SameShapeVectorAllocator for DefaultAllocator where R1: Dim, R2: Dim, - T: Scalar, DefaultAllocator: Allocator + Allocator>, ShapeConstraint: SameNumberOfRows, { diff --git a/src/base/array_storage.rs b/src/base/array_storage.rs index dc4e0df7..3fc88ade 100644 --- a/src/base/array_storage.rs +++ b/src/base/array_storage.rs @@ -12,8 +12,6 @@ use serde::ser::SerializeSeq; use serde::{Deserialize, Deserializer, Serialize, Serializer}; #[cfg(feature = "serde-serialize-no-std")] use std::marker::PhantomData; -#[cfg(feature = "serde-serialize-no-std")] -use std::mem; #[cfg(feature = "abomonation-serialize")] use abomonation::Abomonation; @@ -21,21 +19,37 @@ use abomonation::Abomonation; use crate::base::allocator::Allocator; use crate::base::default_allocator::DefaultAllocator; use crate::base::dimension::{Const, ToTypenum}; -use crate::base::storage::{ - ContiguousStorage, ContiguousStorageMut, Owned, ReshapableStorage, Storage, StorageMut, -}; +use crate::base::storage::{IsContiguous, Owned, RawStorage, RawStorageMut, ReshapableStorage}; use crate::base::Scalar; +use crate::Storage; +use std::mem; /* * - * Static Storage. + * Static RawStorage. * */ /// A array-based statically sized matrix data storage. -#[repr(C)] +#[repr(transparent)] #[derive(Copy, Clone, PartialEq, Eq, Hash)] pub struct ArrayStorage(pub [[T; R]; C]); +impl ArrayStorage { + /// Converts this array storage to a slice. + #[inline] + pub fn as_slice(&self) -> &[T] { + // SAFETY: this is OK because ArrayStorage is contiguous. + unsafe { self.as_slice_unchecked() } + } + + /// Converts this array storage to a mutable slice. + #[inline] + pub fn as_mut_slice(&mut self) -> &mut [T] { + // SAFETY: this is OK because ArrayStorage is contiguous. + unsafe { self.as_mut_slice_unchecked() } + } +} + // TODO: remove this once the stdlib implements Default for arrays. impl Default for ArrayStorage where @@ -54,11 +68,8 @@ impl Debug for ArrayStorage { } } -unsafe impl Storage, Const> +unsafe impl RawStorage, Const> for ArrayStorage -where - T: Scalar, - DefaultAllocator: Allocator, Const, Buffer = Self>, { type RStride = Const<1>; type CStride = Const; @@ -83,6 +94,17 @@ where true } + #[inline] + unsafe fn as_slice_unchecked(&self) -> &[T] { + std::slice::from_raw_parts(self.ptr(), R * C) + } +} + +unsafe impl Storage, Const> + for ArrayStorage +where + DefaultAllocator: Allocator, Const, Buffer = Self>, +{ #[inline] fn into_owned(self) -> Owned, Const> where @@ -96,21 +118,12 @@ where where DefaultAllocator: Allocator, Const>, { - let it = self.as_slice().iter().cloned(); - DefaultAllocator::allocate_from_iterator(self.shape().0, self.shape().1, it) - } - - #[inline] - unsafe fn as_slice_unchecked(&self) -> &[T] { - std::slice::from_raw_parts(self.ptr(), R * C) + self.clone() } } -unsafe impl StorageMut, Const> +unsafe impl RawStorageMut, Const> for ArrayStorage -where - T: Scalar, - DefaultAllocator: Allocator, Const, Buffer = Self>, { #[inline] fn ptr_mut(&mut self) -> *mut T { @@ -123,21 +136,7 @@ where } } -unsafe impl ContiguousStorage, Const> - for ArrayStorage -where - T: Scalar, - DefaultAllocator: Allocator, Const, Buffer = Self>, -{ -} - -unsafe impl ContiguousStorageMut, Const> - for ArrayStorage -where - T: Scalar, - DefaultAllocator: Allocator, Const, Buffer = Self>, -{ -} +unsafe impl IsContiguous for ArrayStorage {} impl ReshapableStorage, Const, Const, Const> for ArrayStorage @@ -160,8 +159,8 @@ where fn reshape_generic(self, _: Const, _: Const) -> Self::Output { unsafe { - let data: [[T; R2]; C2] = std::mem::transmute_copy(&self.0); - std::mem::forget(self.0); + let data: [[T; R2]; C2] = mem::transmute_copy(&self.0); + mem::forget(self.0); ArrayStorage(data) } } @@ -240,19 +239,28 @@ where where V: SeqAccess<'a>, { - let mut out: Self::Value = unsafe { mem::MaybeUninit::uninit().assume_init() }; + let mut out: ArrayStorage, R, C> = + DefaultAllocator::allocate_uninit(Const::, Const::); let mut curr = 0; while let Some(value) = visitor.next_element()? { *out.as_mut_slice() .get_mut(curr) - .ok_or_else(|| V::Error::invalid_length(curr, &self))? = value; + .ok_or_else(|| V::Error::invalid_length(curr, &self))? = + core::mem::MaybeUninit::new(value); curr += 1; } if curr == R * C { - Ok(out) + // Safety: all the elements have been initialized. + unsafe { Ok(, Const>>::assume_init(out)) } } else { + for i in 0..curr { + // Safety: + // - We couldn’t initialize the whole storage. Drop the ones we initialized. + unsafe { std::ptr::drop_in_place(out.as_mut_slice()[i].as_mut_ptr()) }; + } + Err(V::Error::invalid_length(curr, &self)) } } diff --git a/src/base/blas.rs b/src/base/blas.rs index 72b00bda..4d5a5b5d 100644 --- a/src/base/blas.rs +++ b/src/base/blas.rs @@ -1,23 +1,21 @@ -use crate::SimdComplexField; -#[cfg(feature = "std")] -use matrixmultiply; +use crate::{RawStorage, SimdComplexField}; use num::{One, Zero}; use simba::scalar::{ClosedAdd, ClosedMul}; -#[cfg(feature = "std")] -use std::mem; use crate::base::allocator::Allocator; +use crate::base::blas_uninit::{axcpy_uninit, gemm_uninit, gemv_uninit}; use crate::base::constraint::{ AreMultipliable, DimEq, SameNumberOfColumns, SameNumberOfRows, ShapeConstraint, }; use crate::base::dimension::{Const, Dim, Dynamic, U1, U2, U3, U4}; use crate::base::storage::{Storage, StorageMut}; +use crate::base::uninit::Init; use crate::base::{ DVectorSlice, DefaultAllocator, Matrix, Scalar, SquareMatrix, Vector, VectorSlice, }; /// # Dot/scalar product -impl> Matrix +impl> Matrix where T: Scalar + Zero + ClosedAdd + ClosedMul, { @@ -28,7 +26,7 @@ where conjugate: impl Fn(T) -> T, ) -> T where - SB: Storage, + SB: RawStorage, ShapeConstraint: DimEq + DimEq, { assert!( @@ -196,7 +194,7 @@ where #[must_use] pub fn dot(&self, rhs: &Matrix) -> T where - SB: Storage, + SB: RawStorage, ShapeConstraint: DimEq + DimEq, { self.dotx(rhs, |e| e) @@ -226,7 +224,7 @@ where pub fn dotc(&self, rhs: &Matrix) -> T where T: SimdComplexField, - SB: Storage, + SB: RawStorage, ShapeConstraint: DimEq + DimEq, { self.dotx(rhs, T::simd_conjugate) @@ -253,7 +251,7 @@ where #[must_use] pub fn tr_dot(&self, rhs: &Matrix) -> T where - SB: Storage, + SB: RawStorage, ShapeConstraint: DimEq + DimEq, { let (nrows, ncols) = self.shape(); @@ -278,43 +276,6 @@ where } } -#[allow(clippy::too_many_arguments)] -fn array_axcpy( - y: &mut [T], - a: T, - x: &[T], - c: T, - beta: T, - stride1: usize, - stride2: usize, - len: usize, -) where - T: Scalar + Zero + ClosedAdd + ClosedMul, -{ - for i in 0..len { - unsafe { - let y = y.get_unchecked_mut(i * stride1); - *y = a.inlined_clone() - * x.get_unchecked(i * stride2).inlined_clone() - * c.inlined_clone() - + beta.inlined_clone() * y.inlined_clone(); - } - } -} - -fn array_axc(y: &mut [T], a: T, x: &[T], c: T, stride1: usize, stride2: usize, len: usize) -where - T: Scalar + Zero + ClosedAdd + ClosedMul, -{ - for i in 0..len { - unsafe { - *y.get_unchecked_mut(i * stride1) = a.inlined_clone() - * x.get_unchecked(i * stride2).inlined_clone() - * c.inlined_clone(); - } - } -} - /// # BLAS functions impl Vector where @@ -341,23 +302,7 @@ where SB: Storage, ShapeConstraint: DimEq, { - assert_eq!(self.nrows(), x.nrows(), "Axcpy: mismatched vector shapes."); - - let rstride1 = self.strides().0; - let rstride2 = x.strides().0; - - unsafe { - // SAFETY: the conversion to slices is OK because we access the - // elements taking the strides into account. - let y = self.data.as_mut_slice_unchecked(); - let x = x.data.as_slice_unchecked(); - - if !b.is_zero() { - array_axcpy(y, a, x, c, b, rstride1, rstride2, x.len()); - } else { - array_axc(y, a, x, c, rstride1, rstride2, x.len()); - } - } + unsafe { axcpy_uninit(Init, self, a, x, c, b) }; } /// Computes `self = a * x + b * self`. @@ -413,38 +358,8 @@ where SC: Storage, ShapeConstraint: DimEq + AreMultipliable, { - let dim1 = self.nrows(); - let (nrows2, ncols2) = a.shape(); - let dim3 = x.nrows(); - - assert!( - ncols2 == dim3 && dim1 == nrows2, - "Gemv: dimensions mismatch." - ); - - if ncols2 == 0 { - // NOTE: we can't just always multiply by beta - // because we documented the guaranty that `self` is - // never read if `beta` is zero. - if beta.is_zero() { - self.fill(T::zero()); - } else { - *self *= beta; - } - return; - } - - // TODO: avoid bound checks. - let col2 = a.column(0); - let val = unsafe { x.vget_unchecked(0).inlined_clone() }; - self.axcpy(alpha.inlined_clone(), &col2, val, beta); - - for j in 1..ncols2 { - let col2 = a.column(j); - let val = unsafe { x.vget_unchecked(j).inlined_clone() }; - - self.axcpy(alpha.inlined_clone(), &col2, val, T::one()); - } + // Safety: this is safe because we are passing Status == Init. + unsafe { gemv_uninit(Init, self, alpha, a, x, beta) } } #[inline(always)] @@ -504,25 +419,6 @@ where } } - /// Computes `self = alpha * a * x + beta * self`, where `a` is a **symmetric** matrix, `x` a - /// vector, and `alpha, beta` two scalars. DEPRECATED: use `sygemv` instead. - #[inline] - #[deprecated(note = "This is renamed `sygemv` to match the original BLAS terminology.")] - pub fn gemv_symm( - &mut self, - alpha: T, - a: &SquareMatrix, - x: &Vector, - beta: T, - ) where - T: One, - SB: Storage, - SC: Storage, - ShapeConstraint: DimEq + AreMultipliable, - { - self.sygemv(alpha, a, x, beta) - } - /// Computes `self = alpha * a * x + beta * self`, where `a` is a **symmetric** matrix, `x` a /// vector, and `alpha, beta` two scalars. /// @@ -859,122 +755,9 @@ where + SameNumberOfColumns + AreMultipliable, { - let ncols1 = self.ncols(); - - #[cfg(feature = "std")] - { - // We assume large matrices will be Dynamic but small matrices static. - // We could use matrixmultiply for large statically-sized matrices but the performance - // threshold to activate it would be different from SMALL_DIM because our code optimizes - // better for statically-sized matrices. - if R1::is::() - || C1::is::() - || R2::is::() - || C2::is::() - || R3::is::() - || C3::is::() - { - // matrixmultiply can be used only if the std feature is available. - let nrows1 = self.nrows(); - let (nrows2, ncols2) = a.shape(); - let (nrows3, ncols3) = b.shape(); - - // Threshold determined empirically. - const SMALL_DIM: usize = 5; - - if nrows1 > SMALL_DIM - && ncols1 > SMALL_DIM - && nrows2 > SMALL_DIM - && ncols2 > SMALL_DIM - { - assert_eq!( - ncols2, nrows3, - "gemm: dimensions mismatch for multiplication." - ); - assert_eq!( - (nrows1, ncols1), - (nrows2, ncols3), - "gemm: dimensions mismatch for addition." - ); - - // NOTE: this case should never happen because we enter this - // codepath only when ncols2 > SMALL_DIM. Though we keep this - // here just in case if in the future we change the conditions to - // enter this codepath. - if ncols2 == 0 { - // NOTE: we can't just always multiply by beta - // because we documented the guaranty that `self` is - // never read if `beta` is zero. - if beta.is_zero() { - self.fill(T::zero()); - } else { - *self *= beta; - } - return; - } - - if T::is::() { - let (rsa, csa) = a.strides(); - let (rsb, csb) = b.strides(); - let (rsc, csc) = self.strides(); - - unsafe { - matrixmultiply::sgemm( - nrows2, - ncols2, - ncols3, - mem::transmute_copy(&alpha), - a.data.ptr() as *const f32, - rsa as isize, - csa as isize, - b.data.ptr() as *const f32, - rsb as isize, - csb as isize, - mem::transmute_copy(&beta), - self.data.ptr_mut() as *mut f32, - rsc as isize, - csc as isize, - ); - } - return; - } else if T::is::() { - let (rsa, csa) = a.strides(); - let (rsb, csb) = b.strides(); - let (rsc, csc) = self.strides(); - - unsafe { - matrixmultiply::dgemm( - nrows2, - ncols2, - ncols3, - mem::transmute_copy(&alpha), - a.data.ptr() as *const f64, - rsa as isize, - csa as isize, - b.data.ptr() as *const f64, - rsb as isize, - csb as isize, - mem::transmute_copy(&beta), - self.data.ptr_mut() as *mut f64, - rsc as isize, - csc as isize, - ); - } - return; - } - } - } - } - - for j1 in 0..ncols1 { - // TODO: avoid bound checks. - self.column_mut(j1).gemv( - alpha.inlined_clone(), - a, - &b.column(j1), - beta.inlined_clone(), - ); - } + // SAFETY: this is valid because our matrices are initialized and + // we are using status = Init. + unsafe { gemm_uninit(Init, self, alpha, a, b, beta) } } /// Computes `self = alpha * a.transpose() * b + beta * self`, where `a, b, self` are matrices. @@ -1337,9 +1120,8 @@ where ShapeConstraint: DimEq + DimEq + DimEq, DefaultAllocator: Allocator, { - let mut work = unsafe { - crate::unimplemented_or_uninitialized_generic!(self.data.shape().0, Const::<1>) - }; + // TODO: would it be useful to avoid the zero-initialization of the workspace data? + let mut work = Matrix::zeros_generic(self.shape_generic().0, Const::<1>); self.quadform_tr_with_workspace(&mut work, alpha, lhs, mid, beta) } @@ -1432,9 +1214,8 @@ where ShapeConstraint: DimEq + DimEq + AreMultipliable, DefaultAllocator: Allocator, { - let mut work = unsafe { - crate::unimplemented_or_uninitialized_generic!(mid.data.shape().0, Const::<1>) - }; + // TODO: would it be useful to avoid the zero-initialization of the workspace data? + let mut work = Vector::zeros_generic(mid.shape_generic().0, Const::<1>); self.quadform_with_workspace(&mut work, alpha, mid, rhs, beta) } } diff --git a/src/base/blas_uninit.rs b/src/base/blas_uninit.rs new file mode 100644 index 00000000..6f4fde7b --- /dev/null +++ b/src/base/blas_uninit.rs @@ -0,0 +1,323 @@ +/* + * This file implements some BLAS operations in such a way that they work + * even if the first argument (the output parameter) is an uninitialized matrix. + * + * Because doing this makes the code harder to read, we only implemented the operations that we + * know would benefit from this performance-wise, namely, GEMM (which we use for our matrix + * multiplication code). If we identify other operations like that in the future, we could add + * them here. + */ + +#[cfg(feature = "std")] +use matrixmultiply; +use num::{One, Zero}; +use simba::scalar::{ClosedAdd, ClosedMul}; +#[cfg(feature = "std")] +use std::mem; + +use crate::base::constraint::{ + AreMultipliable, DimEq, SameNumberOfColumns, SameNumberOfRows, ShapeConstraint, +}; +use crate::base::dimension::{Dim, Dynamic, U1}; +use crate::base::storage::{RawStorage, RawStorageMut}; +use crate::base::uninit::InitStatus; +use crate::base::{Matrix, Scalar, Vector}; +use std::any::TypeId; + +// # Safety +// The content of `y` must only contain values for which +// `Status::assume_init_mut` is sound. +#[allow(clippy::too_many_arguments)] +unsafe fn array_axcpy( + _: Status, + y: &mut [Status::Value], + a: T, + x: &[T], + c: T, + beta: T, + stride1: usize, + stride2: usize, + len: usize, +) where + Status: InitStatus, + T: Scalar + Zero + ClosedAdd + ClosedMul, +{ + for i in 0..len { + let y = Status::assume_init_mut(y.get_unchecked_mut(i * stride1)); + *y = a.inlined_clone() * x.get_unchecked(i * stride2).inlined_clone() * c.inlined_clone() + + beta.inlined_clone() * y.inlined_clone(); + } +} + +fn array_axc( + _: Status, + y: &mut [Status::Value], + a: T, + x: &[T], + c: T, + stride1: usize, + stride2: usize, + len: usize, +) where + Status: InitStatus, + T: Scalar + Zero + ClosedAdd + ClosedMul, +{ + for i in 0..len { + unsafe { + Status::init( + y.get_unchecked_mut(i * stride1), + a.inlined_clone() + * x.get_unchecked(i * stride2).inlined_clone() + * c.inlined_clone(), + ); + } + } +} + +/// Computes `y = a * x * c + b * y`. +/// +/// If `b` is zero, `y` is never read from and may be uninitialized. +/// +/// # Safety +/// This is UB if b != 0 and any component of `y` is uninitialized. +#[inline(always)] +#[allow(clippy::many_single_char_names)] +pub unsafe fn axcpy_uninit( + status: Status, + y: &mut Vector, + a: T, + x: &Vector, + c: T, + b: T, +) where + T: Scalar + Zero + ClosedAdd + ClosedMul, + SA: RawStorageMut, + SB: RawStorage, + ShapeConstraint: DimEq, + Status: InitStatus, +{ + assert_eq!(y.nrows(), x.nrows(), "Axcpy: mismatched vector shapes."); + + 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(); + + if !b.is_zero() { + array_axcpy(status, y, a, x, c, b, rstride1, rstride2, x.len()); + } else { + array_axc(status, y, a, x, c, rstride1, rstride2, x.len()); + } +} + +/// Computes `y = alpha * a * x + beta * y`, where `a` is a matrix, `x` a vector, and +/// `alpha, beta` two scalars. +/// +/// If `beta` is zero, `y` is never read from and may be uninitialized. +/// +/// # Safety +/// This is UB if beta != 0 and any component of `y` is uninitialized. +#[inline(always)] +pub unsafe fn gemv_uninit( + status: Status, + y: &mut Vector, + alpha: T, + a: &Matrix, + x: &Vector, + beta: T, +) where + Status: InitStatus, + T: Scalar + Zero + One + ClosedAdd + ClosedMul, + SA: RawStorageMut, + SB: RawStorage, + SC: RawStorage, + ShapeConstraint: DimEq + AreMultipliable, +{ + let dim1 = y.nrows(); + let (nrows2, ncols2) = a.shape(); + let dim3 = x.nrows(); + + assert!( + ncols2 == dim3 && dim1 == nrows2, + "Gemv: dimensions mismatch." + ); + + if ncols2 == 0 { + if beta.is_zero() { + y.apply(|e| Status::init(e, T::zero())); + } else { + // SAFETY: this is UB if y is uninitialized. + y.apply(|e| *Status::assume_init_mut(e) *= beta.inlined_clone()); + } + return; + } + + // TODO: avoid bound checks. + let col2 = a.column(0); + let val = x.vget_unchecked(0).inlined_clone(); + + // SAFETY: this is the call that makes this method unsafe: it is UB if Status = Uninit and beta != 0. + axcpy_uninit(status, y, alpha.inlined_clone(), &col2, val, beta); + + for j in 1..ncols2 { + let col2 = a.column(j); + let val = x.vget_unchecked(j).inlined_clone(); + + // SAFETY: safe because y was initialized above. + axcpy_uninit(status, y, alpha.inlined_clone(), &col2, val, T::one()); + } +} + +/// Computes `y = alpha * a * b + beta * y`, where `a, b, y` are matrices. +/// `alpha` and `beta` are scalar. +/// +/// If `beta` is zero, `y` is never read from and may be uninitialized. +/// +/// # Safety +/// This is UB if beta != 0 and any component of `y` is uninitialized. +#[inline(always)] +pub unsafe fn gemm_uninit< + Status, + T, + R1: Dim, + C1: Dim, + R2: Dim, + C2: Dim, + R3: Dim, + C3: Dim, + SA, + SB, + SC, +>( + status: Status, + y: &mut Matrix, + alpha: T, + a: &Matrix, + b: &Matrix, + beta: T, +) where + Status: InitStatus, + T: Scalar + Zero + One + ClosedAdd + ClosedMul, + SA: RawStorageMut, + SB: RawStorage, + SC: RawStorage, + ShapeConstraint: + SameNumberOfRows + SameNumberOfColumns + AreMultipliable, +{ + let ncols1 = y.ncols(); + + #[cfg(feature = "std")] + { + // We assume large matrices will be Dynamic but small matrices static. + // We could use matrixmultiply for large statically-sized matrices but the performance + // threshold to activate it would be different from SMALL_DIM because our code optimizes + // better for statically-sized matrices. + if R1::is::() + || C1::is::() + || R2::is::() + || C2::is::() + || R3::is::() + || C3::is::() + { + // matrixmultiply can be used only if the std feature is available. + let nrows1 = y.nrows(); + let (nrows2, ncols2) = a.shape(); + let (nrows3, ncols3) = b.shape(); + + // Threshold determined empirically. + const SMALL_DIM: usize = 5; + + if nrows1 > SMALL_DIM && ncols1 > SMALL_DIM && nrows2 > SMALL_DIM && ncols2 > SMALL_DIM + { + assert_eq!( + ncols2, nrows3, + "gemm: dimensions mismatch for multiplication." + ); + assert_eq!( + (nrows1, ncols1), + (nrows2, ncols3), + "gemm: dimensions mismatch for addition." + ); + + // NOTE: this case should never happen because we enter this + // codepath only when ncols2 > SMALL_DIM. Though we keep this + // here just in case if in the future we change the conditions to + // enter this codepath. + if ncols2 == 0 { + // NOTE: we can't just always multiply by beta + // because we documented the guaranty that `self` is + // never read if `beta` is zero. + if beta.is_zero() { + y.apply(|e| Status::init(e, T::zero())); + } else { + // SAFETY: this is UB if Status = Uninit + y.apply(|e| *Status::assume_init_mut(e) *= beta.inlined_clone()); + } + return; + } + + if TypeId::of::() == TypeId::of::() { + let (rsa, csa) = a.strides(); + let (rsb, csb) = b.strides(); + let (rsc, csc) = y.strides(); + + matrixmultiply::sgemm( + nrows2, + ncols2, + ncols3, + mem::transmute_copy(&alpha), + a.data.ptr() as *const f32, + rsa as isize, + csa as isize, + b.data.ptr() as *const f32, + rsb as isize, + csb as isize, + mem::transmute_copy(&beta), + y.data.ptr_mut() as *mut f32, + rsc as isize, + csc as isize, + ); + return; + } else if TypeId::of::() == TypeId::of::() { + let (rsa, csa) = a.strides(); + let (rsb, csb) = b.strides(); + let (rsc, csc) = y.strides(); + + matrixmultiply::dgemm( + nrows2, + ncols2, + ncols3, + mem::transmute_copy(&alpha), + a.data.ptr() as *const f64, + rsa as isize, + csa as isize, + b.data.ptr() as *const f64, + rsb as isize, + csb as isize, + mem::transmute_copy(&beta), + y.data.ptr_mut() as *mut f64, + rsc as isize, + csc as isize, + ); + return; + } + } + } + } + + for j1 in 0..ncols1 { + // TODO: avoid bound checks. + // SAFETY: this is UB if Status = Uninit && beta != 0 + gemv_uninit( + status, + &mut y.column_mut(j1), + alpha.inlined_clone(), + a, + &b.column(j1), + beta.inlined_clone(), + ); + } +} diff --git a/src/base/construction.rs b/src/base/construction.rs index cde4f924..3deb66c2 100644 --- a/src/base/construction.rs +++ b/src/base/construction.rs @@ -14,33 +14,32 @@ use rand::{ }; use std::iter; -use std::mem; use typenum::{self, Cmp, Greater}; use simba::scalar::{ClosedAdd, ClosedMul}; use crate::base::allocator::Allocator; use crate::base::dimension::{Dim, DimName, Dynamic, ToTypenum}; -use crate::base::storage::Storage; +use crate::base::storage::RawStorage; use crate::base::{ ArrayStorage, Const, DefaultAllocator, Matrix, OMatrix, OVector, Scalar, Unit, Vector, }; +use crate::UninitMatrix; +use std::mem::MaybeUninit; -/// When "no_unsound_assume_init" is enabled, expands to `unimplemented!()` instead of `new_uninitialized_generic().assume_init()`. -/// Intended as a placeholder, each callsite should be refactored to use uninitialized memory soundly -#[macro_export] -macro_rules! unimplemented_or_uninitialized_generic { - ($nrows:expr, $ncols:expr) => {{ - #[cfg(feature="no_unsound_assume_init")] { - // Some of the call sites need the number of rows and columns from this to infer a type, so - // uninitialized memory is used to infer the type, as `T: Zero` isn't available at all callsites. - // This may technically still be UB even though the assume_init is dead code, but all callsites should be fixed before #556 is closed. - let typeinference_helper = crate::base::Matrix::new_uninitialized_generic($nrows, $ncols); - unimplemented!(); - typeinference_helper.assume_init() +impl UninitMatrix +where + DefaultAllocator: Allocator, +{ + /// Builds a matrix with uninitialized elements of type `MaybeUninit`. + #[inline(always)] + pub fn uninit(nrows: R, ncols: C) -> Self { + // SAFETY: this is OK because the dimension automatically match the storage + // because we are building an owned storage. + unsafe { + Self::from_data_statically_unchecked(DefaultAllocator::allocate_uninit(nrows, ncols)) } - #[cfg(not(feature="no_unsound_assume_init"))] { crate::base::Matrix::new_uninitialized_generic($nrows, $ncols).assume_init() } - }} + } } /// # Generic constructors @@ -53,16 +52,6 @@ impl OMatrix where DefaultAllocator: Allocator, { - /// Creates a new uninitialized matrix. - /// - /// # Safety - /// If the matrix has a compile-time dimension, this panics - /// if `nrows != R::to_usize()` or `ncols != C::to_usize()`. - #[inline] - pub unsafe fn new_uninitialized_generic(nrows: R, ncols: C) -> mem::MaybeUninit { - Self::from_uninitialized_data(DefaultAllocator::allocate_uninitialized(nrows, ncols)) - } - /// Creates a matrix with all its elements set to `elem`. #[inline] pub fn from_element_generic(nrows: R, ncols: C, elem: T) -> Self { @@ -109,16 +98,20 @@ where "Matrix init. error: the slice did not contain the right number of elements." ); - let mut res = unsafe { crate::unimplemented_or_uninitialized_generic!(nrows, ncols) }; + let mut res = Matrix::uninit(nrows, ncols); let mut iter = slice.iter(); - for i in 0..nrows.value() { - for j in 0..ncols.value() { - unsafe { *res.get_unchecked_mut((i, j)) = iter.next().unwrap().inlined_clone() } + unsafe { + for i in 0..nrows.value() { + for j in 0..ncols.value() { + *res.get_unchecked_mut((i, j)) = + MaybeUninit::new(iter.next().unwrap().inlined_clone()) + } } - } - res + // SAFETY: the result has been fully initialized above. + res.assume_init() + } } /// Creates a matrix with its elements filled with the components provided by a slice. The @@ -135,15 +128,18 @@ where where F: FnMut(usize, usize) -> T, { - let mut res: Self = unsafe { crate::unimplemented_or_uninitialized_generic!(nrows, ncols) }; + let mut res = Matrix::uninit(nrows, ncols); - for j in 0..ncols.value() { - for i in 0..nrows.value() { - unsafe { *res.get_unchecked_mut((i, j)) = f(i, j) } + unsafe { + for j in 0..ncols.value() { + for i in 0..nrows.value() { + *res.get_unchecked_mut((i, j)) = MaybeUninit::new(f(i, j)); + } } - } - res + // SAFETY: the result has been fully initialized above. + res.assume_init() + } } /// Creates a new identity matrix. @@ -217,7 +213,7 @@ where #[inline] pub fn from_rows(rows: &[Matrix, C, SB>]) -> Self where - SB: Storage, C>, + SB: RawStorage, C>, { assert!(!rows.is_empty(), "At least one row must be given."); let nrows = R::try_to_usize().unwrap_or_else(|| rows.len()); @@ -259,7 +255,7 @@ where #[inline] pub fn from_columns(columns: &[Vector]) -> Self where - SB: Storage, + SB: RawStorage, { assert!(!columns.is_empty(), "At least one column must be given."); let ncols = C::try_to_usize().unwrap_or_else(|| columns.len()); @@ -353,11 +349,11 @@ where /// dm[(2, 0)] == 0.0 && dm[(2, 1)] == 0.0 && dm[(2, 2)] == 3.0); /// ``` #[inline] - pub fn from_diagonal>(diag: &Vector) -> Self + pub fn from_diagonal>(diag: &Vector) -> Self where T: Zero, { - let (dim, _) = diag.data.shape(); + let (dim, _) = diag.shape_generic(); let mut res = Self::zeros_generic(dim, dim); for i in 0..diag.len() { @@ -377,12 +373,6 @@ where */ macro_rules! impl_constructors( ($($Dims: ty),*; $(=> $DimIdent: ident: $DimBound: ident),*; $($gargs: expr),*; $($args: ident),*) => { - /// Creates a new uninitialized matrix or vector. - #[inline] - pub unsafe fn new_uninitialized($($args: usize),*) -> mem::MaybeUninit { - Self::new_uninitialized_generic($($gargs),*) - } - /// Creates a matrix or vector with all its elements set to `elem`. /// /// # Example diff --git a/src/base/conversion.rs b/src/base/conversion.rs index 8ede11ca..ec7fd936 100644 --- a/src/base/conversion.rs +++ b/src/base/conversion.rs @@ -14,7 +14,7 @@ use crate::base::dimension::{ Const, Dim, DimName, U1, U10, U11, U12, U13, U14, U15, U16, U2, U3, U4, U5, U6, U7, U8, U9, }; use crate::base::iter::{MatrixIter, MatrixIterMut}; -use crate::base::storage::{ContiguousStorage, ContiguousStorageMut, Storage, StorageMut}; +use crate::base::storage::{IsContiguous, RawStorage, RawStorageMut}; use crate::base::{ ArrayStorage, DVectorSlice, DVectorSliceMut, DefaultAllocator, Matrix, MatrixSlice, MatrixSliceMut, OMatrix, Scalar, @@ -24,6 +24,7 @@ use crate::base::{DVector, VecStorage}; use crate::base::{SliceStorage, SliceStorageMut}; use crate::constraint::DimEq; use crate::{IsNotStaticOne, RowSVector, SMatrix, SVector}; +use std::mem::MaybeUninit; // TODO: too bad this won't work for slice conversions. impl SubsetOf> for OMatrix @@ -43,18 +44,20 @@ where let (nrows, ncols) = self.shape(); let nrows2 = R2::from_usize(nrows); let ncols2 = C2::from_usize(ncols); + let mut res = Matrix::uninit(nrows2, ncols2); - let mut res: OMatrix = - unsafe { crate::unimplemented_or_uninitialized_generic!(nrows2, ncols2) }; for i in 0..nrows { for j in 0..ncols { + // Safety: all indices are in range. unsafe { - *res.get_unchecked_mut((i, j)) = T2::from_subset(self.get_unchecked((i, j))) + *res.get_unchecked_mut((i, j)) = + MaybeUninit::new(T2::from_subset(self.get_unchecked((i, j)))); } } } - res + // Safety: res is now fully initialized. + unsafe { res.assume_init() } } #[inline] @@ -67,21 +70,25 @@ where let (nrows2, ncols2) = m.shape(); let nrows = R1::from_usize(nrows2); let ncols = C1::from_usize(ncols2); + let mut res = Matrix::uninit(nrows, ncols); - let mut res: Self = unsafe { crate::unimplemented_or_uninitialized_generic!(nrows, ncols) }; for i in 0..nrows2 { for j in 0..ncols2 { + // Safety: all indices are in range. unsafe { - *res.get_unchecked_mut((i, j)) = m.get_unchecked((i, j)).to_subset_unchecked() + *res.get_unchecked_mut((i, j)) = + MaybeUninit::new(m.get_unchecked((i, j)).to_subset_unchecked()) } } } - res + unsafe { res.assume_init() } } } -impl<'a, T: Scalar, R: Dim, C: Dim, S: Storage> IntoIterator for &'a Matrix { +impl<'a, T: Scalar, R: Dim, C: Dim, S: RawStorage> IntoIterator + for &'a Matrix +{ type Item = &'a T; type IntoIter = MatrixIter<'a, T, R, C, S>; @@ -91,7 +98,7 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: Storage> IntoIterator for &'a Ma } } -impl<'a, T: Scalar, R: Dim, C: Dim, S: StorageMut> IntoIterator +impl<'a, T: Scalar, R: Dim, C: Dim, S: RawStorageMut> IntoIterator for &'a mut Matrix { type Item = &'a mut T; @@ -142,9 +149,10 @@ macro_rules! impl_from_into_asref_1D( ($(($NRows: ident, $NCols: ident) => $SZ: expr);* $(;)*) => {$( impl AsRef<[T; $SZ]> for Matrix where T: Scalar, - S: ContiguousStorage { + S: RawStorage + IsContiguous { #[inline] fn as_ref(&self) -> &[T; $SZ] { + // Safety: this is OK thanks to the IsContiguous trait. unsafe { &*(self.data.ptr() as *const [T; $SZ]) } @@ -153,9 +161,10 @@ macro_rules! impl_from_into_asref_1D( impl AsMut<[T; $SZ]> for Matrix where T: Scalar, - S: ContiguousStorageMut { + S: RawStorageMut + IsContiguous { #[inline] fn as_mut(&mut self) -> &mut [T; $SZ] { + // Safety: this is OK thanks to the IsContiguous trait. unsafe { &mut *(self.data.ptr_mut() as *mut [T; $SZ]) } @@ -201,9 +210,10 @@ macro_rules! impl_from_into_asref_borrow_2D( $Ref:ident.$ref:ident(), $Mut:ident.$mut:ident() ) => { impl $Ref<[[T; $SZRows]; $SZCols]> for Matrix - where S: ContiguousStorage { + where S: RawStorage + IsContiguous { #[inline] fn $ref(&self) -> &[[T; $SZRows]; $SZCols] { + // Safety: OK thanks to the IsContiguous trait. unsafe { &*(self.data.ptr() as *const [[T; $SZRows]; $SZCols]) } @@ -211,9 +221,10 @@ macro_rules! impl_from_into_asref_borrow_2D( } impl $Mut<[[T; $SZRows]; $SZCols]> for Matrix - where S: ContiguousStorageMut { + where S: RawStorageMut + IsContiguous { #[inline] fn $mut(&mut self) -> &mut [[T; $SZRows]; $SZCols] { + // Safety: OK thanks to the IsContiguous trait. unsafe { &mut *(self.data.ptr_mut() as *mut [[T; $SZRows]; $SZCols]) } @@ -333,14 +344,14 @@ where CSlice: Dim, RStride: Dim, CStride: Dim, - S: Storage, + S: RawStorage, ShapeConstraint: DimEq + DimEq + DimEq + DimEq, { fn from(m: &'a Matrix) -> Self { - let (row, col) = m.data.shape(); + let (row, col) = m.shape_generic(); let row_slice = RSlice::from_usize(row.value()); let col_slice = CSlice::from_usize(col.value()); @@ -370,14 +381,14 @@ where CSlice: Dim, RStride: Dim, CStride: Dim, - S: Storage, + S: RawStorage, ShapeConstraint: DimEq + DimEq + DimEq + DimEq, { fn from(m: &'a mut Matrix) -> Self { - let (row, col) = m.data.shape(); + let (row, col) = m.shape_generic(); let row_slice = RSlice::from_usize(row.value()); let col_slice = CSlice::from_usize(col.value()); @@ -407,14 +418,14 @@ where CSlice: Dim, RStride: Dim, CStride: Dim, - S: StorageMut, + S: RawStorageMut, ShapeConstraint: DimEq + DimEq + DimEq + DimEq, { fn from(m: &'a mut Matrix) -> Self { - let (row, col) = m.data.shape(); + let (row, col) = m.shape_generic(); let row_slice = RSlice::from_usize(row.value()); let col_slice = CSlice::from_usize(col.value()); @@ -442,7 +453,7 @@ impl<'a, T: Scalar> From> for DVector { } } -impl<'a, T: Scalar + Copy, R: Dim, C: Dim, S: ContiguousStorage> +impl<'a, T: Scalar + Copy, R: Dim, C: Dim, S: RawStorage + IsContiguous> From<&'a Matrix> for &'a [T] { #[inline] @@ -451,7 +462,7 @@ impl<'a, T: Scalar + Copy, R: Dim, C: Dim, S: ContiguousStorage> } } -impl<'a, T: Scalar + Copy, R: Dim, C: Dim, S: ContiguousStorageMut> +impl<'a, T: Scalar + Copy, R: Dim, C: Dim, S: RawStorageMut + IsContiguous> From<&'a mut Matrix> for &'a mut [T] { #[inline] @@ -495,7 +506,7 @@ where { #[inline] fn from(arr: [OMatrix; 2]) -> Self { - let (nrows, ncols) = arr[0].data.shape(); + let (nrows, ncols) = arr[0].shape_generic(); Self::from_fn_generic(nrows, ncols, |i, j| { [ @@ -516,7 +527,7 @@ where { #[inline] fn from(arr: [OMatrix; 4]) -> Self { - let (nrows, ncols) = arr[0].data.shape(); + let (nrows, ncols) = arr[0].shape_generic(); Self::from_fn_generic(nrows, ncols, |i, j| { [ @@ -539,7 +550,7 @@ where { #[inline] fn from(arr: [OMatrix; 8]) -> Self { - let (nrows, ncols) = arr[0].data.shape(); + let (nrows, ncols) = arr[0].shape_generic(); Self::from_fn_generic(nrows, ncols, |i, j| { [ @@ -565,7 +576,7 @@ where DefaultAllocator: Allocator + Allocator, { fn from(arr: [OMatrix; 16]) -> Self { - let (nrows, ncols) = arr[0].data.shape(); + let (nrows, ncols) = arr[0].shape_generic(); Self::from_fn_generic(nrows, ncols, |i, j| { [ diff --git a/src/base/coordinates.rs b/src/base/coordinates.rs index be05d3e5..db66811d 100644 --- a/src/base/coordinates.rs +++ b/src/base/coordinates.rs @@ -7,7 +7,7 @@ use std::ops::{Deref, DerefMut}; use crate::base::dimension::{U1, U2, U3, U4, U5, U6}; -use crate::base::storage::{ContiguousStorage, ContiguousStorageMut}; +use crate::base::storage::{IsContiguous, RawStorage, RawStorageMut}; use crate::base::{Matrix, Scalar}; /* @@ -32,19 +32,21 @@ macro_rules! coords_impl( macro_rules! deref_impl( ($R: ty, $C: ty; $Target: ident) => { impl Deref for Matrix - where S: ContiguousStorage { + where S: RawStorage + IsContiguous { type Target = $Target; #[inline] fn deref(&self) -> &Self::Target { + // Safety: this is OK because of the IsContiguous trait. unsafe { &*(self.data.ptr() as *const Self::Target) } } } impl DerefMut for Matrix - where S: ContiguousStorageMut { + where S: RawStorageMut + IsContiguous { #[inline] fn deref_mut(&mut self) -> &mut Self::Target { + // Safety: this is OK because of the IsContiguous trait. unsafe { &mut *(self.data.ptr_mut() as *mut Self::Target) } } } diff --git a/src/base/default_allocator.rs b/src/base/default_allocator.rs index b053c829..b676b5e3 100644 --- a/src/base/default_allocator.rs +++ b/src/base/default_allocator.rs @@ -4,7 +4,6 @@ //! heap-allocated buffers for matrices with at least one dimension unknown at compile-time. use std::cmp; -use std::mem; use std::ptr; #[cfg(all(feature = "alloc", not(feature = "std")))] @@ -16,10 +15,11 @@ use crate::base::array_storage::ArrayStorage; #[cfg(any(feature = "alloc", feature = "std"))] use crate::base::dimension::Dynamic; use crate::base::dimension::{Dim, DimName}; -use crate::base::storage::{ContiguousStorageMut, Storage, StorageMut}; +use crate::base::storage::{RawStorage, RawStorageMut}; #[cfg(any(feature = "std", feature = "alloc"))] use crate::base::vec_storage::VecStorage; use crate::base::Scalar; +use std::mem::{ManuallyDrop, MaybeUninit}; /* * @@ -36,10 +36,23 @@ impl Allocator, Const> for DefaultAllocator { type Buffer = ArrayStorage; + type BufferUninit = ArrayStorage, R, C>; - #[inline] - unsafe fn allocate_uninitialized(_: Const, _: Const) -> mem::MaybeUninit { - mem::MaybeUninit::::uninit() + #[inline(always)] + fn allocate_uninit(_: Const, _: Const) -> ArrayStorage, R, C> { + // SAFETY: An uninitialized `[MaybeUninit<_>; _]` is valid. + let array: [[MaybeUninit; R]; C] = unsafe { MaybeUninit::uninit().assume_init() }; + ArrayStorage(array) + } + + #[inline(always)] + unsafe fn assume_init(uninit: ArrayStorage, R, C>) -> ArrayStorage { + // Safety: + // * The caller guarantees that all elements of the array are initialized + // * `MaybeUninit` and T are guaranteed to have the same layout + // * `MaybeUninit` does not drop, so there are no double-frees + // And thus the conversion is safe + ArrayStorage((&uninit as *const _ as *const [_; C]).read()) } #[inline] @@ -48,14 +61,13 @@ impl Allocator, Const> ncols: Const, iter: I, ) -> Self::Buffer { - #[cfg(feature = "no_unsound_assume_init")] - let mut res: Self::Buffer = unimplemented!(); - #[cfg(not(feature = "no_unsound_assume_init"))] - let mut res = unsafe { Self::allocate_uninitialized(nrows, ncols).assume_init() }; + let mut res = Self::allocate_uninit(nrows, ncols); let mut count = 0; - for (res, e) in res.as_mut_slice().iter_mut().zip(iter.into_iter()) { - *res = e; + // Safety: conversion to a slice is OK because the Buffer is known to be contiguous. + let res_slice = unsafe { res.as_mut_slice_unchecked() }; + for (res, e) in res_slice.iter_mut().zip(iter.into_iter()) { + *res = MaybeUninit::new(e); count += 1; } @@ -64,7 +76,9 @@ impl Allocator, Const> "Matrix init. from iterator: iterator not long enough." ); - res + // Safety: the assertion above made sure that the iterator + // yielded enough elements to initialize our matrix. + unsafe { , Const>>::assume_init(res) } } } @@ -73,15 +87,32 @@ impl Allocator, Const> #[cfg(any(feature = "std", feature = "alloc"))] impl Allocator for DefaultAllocator { type Buffer = VecStorage; + type BufferUninit = VecStorage, Dynamic, C>; #[inline] - unsafe fn allocate_uninitialized(nrows: Dynamic, ncols: C) -> mem::MaybeUninit { - let mut res = Vec::new(); + fn allocate_uninit(nrows: Dynamic, ncols: C) -> VecStorage, Dynamic, C> { + let mut data = Vec::new(); let length = nrows.value() * ncols.value(); - res.reserve_exact(length); - res.set_len(length); + data.reserve_exact(length); + data.resize_with(length, MaybeUninit::uninit); + VecStorage::new(nrows, ncols, data) + } - mem::MaybeUninit::new(VecStorage::new(nrows, ncols, res)) + #[inline] + unsafe fn assume_init( + uninit: VecStorage, Dynamic, C>, + ) -> VecStorage { + // Avoids a double-drop. + let (nrows, ncols) = uninit.shape(); + let vec: Vec<_> = uninit.into(); + let mut md = ManuallyDrop::new(vec); + + // Safety: + // - MaybeUninit has the same alignment and layout as T. + // - The length and capacity come from a valid vector. + let new_data = Vec::from_raw_parts(md.as_mut_ptr() as *mut _, md.len(), md.capacity()); + + VecStorage::new(nrows, ncols, new_data) } #[inline] @@ -103,15 +134,33 @@ impl Allocator for DefaultAllocator { #[cfg(any(feature = "std", feature = "alloc"))] impl Allocator for DefaultAllocator { type Buffer = VecStorage; + type BufferUninit = VecStorage, R, Dynamic>; #[inline] - unsafe fn allocate_uninitialized(nrows: R, ncols: Dynamic) -> mem::MaybeUninit { - let mut res = Vec::new(); + fn allocate_uninit(nrows: R, ncols: Dynamic) -> VecStorage, R, Dynamic> { + let mut data = Vec::new(); let length = nrows.value() * ncols.value(); - res.reserve_exact(length); - res.set_len(length); + data.reserve_exact(length); + data.resize_with(length, MaybeUninit::uninit); - mem::MaybeUninit::new(VecStorage::new(nrows, ncols, res)) + VecStorage::new(nrows, ncols, data) + } + + #[inline] + unsafe fn assume_init( + uninit: VecStorage, R, Dynamic>, + ) -> VecStorage { + // Avoids a double-drop. + let (nrows, ncols) = uninit.shape(); + let vec: Vec<_> = uninit.into(); + let mut md = ManuallyDrop::new(vec); + + // Safety: + // - MaybeUninit has the same alignment and layout as T. + // - The length and capacity come from a valid vector. + let new_data = Vec::from_raw_parts(md.as_mut_ptr() as *mut _, md.len(), md.capacity()); + + VecStorage::new(nrows, ncols, new_data) } #[inline] @@ -146,20 +195,21 @@ where unsafe fn reallocate_copy( rto: Const, cto: Const, - buf: >::Buffer, - ) -> ArrayStorage { - #[cfg(feature = "no_unsound_assume_init")] - let mut res: ArrayStorage = unimplemented!(); - #[cfg(not(feature = "no_unsound_assume_init"))] - let mut res = - , Const>>::allocate_uninitialized(rto, cto) - .assume_init(); + mut buf: >::Buffer, + ) -> ArrayStorage, RTO, CTO> { + let mut res = , Const>>::allocate_uninit(rto, cto); let (rfrom, cfrom) = buf.shape(); let len_from = rfrom.value() * cfrom.value(); let len_to = rto.value() * cto.value(); - ptr::copy_nonoverlapping(buf.ptr(), res.ptr_mut(), cmp::min(len_from, len_to)); + let len_copied = cmp::min(len_from, len_to); + ptr::copy_nonoverlapping(buf.ptr(), res.ptr_mut() as *mut T, len_copied); + + // Safety: + // - We don’t care about dropping elements because the caller is responsible for dropping things. + // - We forget `buf` so that we don’t drop the other elements. + std::mem::forget(buf); res } @@ -176,19 +226,21 @@ where unsafe fn reallocate_copy( rto: Dynamic, cto: CTo, - buf: ArrayStorage, - ) -> VecStorage { - #[cfg(feature = "no_unsound_assume_init")] - let mut res: VecStorage = unimplemented!(); - #[cfg(not(feature = "no_unsound_assume_init"))] - let mut res = - >::allocate_uninitialized(rto, cto).assume_init(); + mut buf: ArrayStorage, + ) -> VecStorage, Dynamic, CTo> { + let mut res = >::allocate_uninit(rto, cto); let (rfrom, cfrom) = buf.shape(); let len_from = rfrom.value() * cfrom.value(); let len_to = rto.value() * cto.value(); - ptr::copy_nonoverlapping(buf.ptr(), res.ptr_mut(), cmp::min(len_from, len_to)); + let len_copied = cmp::min(len_from, len_to); + ptr::copy_nonoverlapping(buf.ptr(), res.ptr_mut() as *mut T, len_copied); + + // Safety: + // - We don’t care about dropping elements because the caller is responsible for dropping things. + // - We forget `buf` so that we don’t drop the other elements. + std::mem::forget(buf); res } @@ -205,19 +257,21 @@ where unsafe fn reallocate_copy( rto: RTo, cto: Dynamic, - buf: ArrayStorage, - ) -> VecStorage { - #[cfg(feature = "no_unsound_assume_init")] - let mut res: VecStorage = unimplemented!(); - #[cfg(not(feature = "no_unsound_assume_init"))] - let mut res = - >::allocate_uninitialized(rto, cto).assume_init(); + mut buf: ArrayStorage, + ) -> VecStorage, RTo, Dynamic> { + let mut res = >::allocate_uninit(rto, cto); let (rfrom, cfrom) = buf.shape(); let len_from = rfrom.value() * cfrom.value(); let len_to = rto.value() * cto.value(); - ptr::copy_nonoverlapping(buf.ptr(), res.ptr_mut(), cmp::min(len_from, len_to)); + let len_copied = cmp::min(len_from, len_to); + ptr::copy_nonoverlapping(buf.ptr(), res.ptr_mut() as *mut T, len_copied); + + // Safety: + // - We don’t care about dropping elements because the caller is responsible for dropping things. + // - We forget `buf` so that we don’t drop the other elements. + std::mem::forget(buf); res } @@ -233,7 +287,7 @@ impl Reallocator, - ) -> VecStorage { + ) -> VecStorage, Dynamic, CTo> { let new_buf = buf.resize(rto.value() * cto.value()); VecStorage::new(rto, cto, new_buf) } @@ -248,7 +302,7 @@ impl Reallocator, - ) -> VecStorage { + ) -> VecStorage, RTo, Dynamic> { let new_buf = buf.resize(rto.value() * cto.value()); VecStorage::new(rto, cto, new_buf) } @@ -263,7 +317,7 @@ impl Reallocator, - ) -> VecStorage { + ) -> VecStorage, Dynamic, CTo> { let new_buf = buf.resize(rto.value() * cto.value()); VecStorage::new(rto, cto, new_buf) } @@ -278,7 +332,7 @@ impl Reallocator, - ) -> VecStorage { + ) -> VecStorage, RTo, Dynamic> { let new_buf = buf.resize(rto.value() * cto.value()); VecStorage::new(rto, cto, new_buf) } diff --git a/src/base/edition.rs b/src/base/edition.rs index f403f9d3..bca017c4 100644 --- a/src/base/edition.rs +++ b/src/base/edition.rs @@ -2,8 +2,6 @@ use num::{One, Zero}; use std::cmp; #[cfg(any(feature = "std", feature = "alloc"))] use std::iter::ExactSizeIterator; -#[cfg(any(feature = "std", feature = "alloc"))] -use std::mem; use std::ptr; use crate::base::allocator::{Allocator, Reallocator}; @@ -11,8 +9,10 @@ use crate::base::constraint::{DimEq, SameNumberOfColumns, SameNumberOfRows, Shap #[cfg(any(feature = "std", feature = "alloc"))] use crate::base::dimension::Dynamic; use crate::base::dimension::{Const, Dim, DimAdd, DimDiff, DimMin, DimMinimum, DimSub, DimSum, U1}; -use crate::base::storage::{ContiguousStorageMut, ReshapableStorage, Storage, StorageMut}; +use crate::base::storage::{RawStorage, RawStorageMut, ReshapableStorage}; use crate::base::{DefaultAllocator, Matrix, OMatrix, RowVector, Scalar, Vector}; +use crate::{Storage, UninitMatrix}; +use std::mem::MaybeUninit; /// # Rows and columns extraction impl> Matrix { @@ -52,10 +52,8 @@ impl> Matrix { DefaultAllocator: Allocator, { let irows = irows.into_iter(); - let ncols = self.data.shape().1; - let mut res = unsafe { - crate::unimplemented_or_uninitialized_generic!(Dynamic::new(irows.len()), ncols) - }; + let ncols = self.shape_generic().1; + let mut res = Matrix::uninit(Dynamic::new(irows.len()), ncols); // First, check that all the indices from irows are valid. // This will allow us to use unchecked access in the inner loop. @@ -69,14 +67,16 @@ impl> Matrix { let src = self.column(j); for (destination, source) in irows.clone().enumerate() { + // Safety: all indices are in range. unsafe { *res.vget_unchecked_mut(destination) = - src.vget_unchecked(*source).inlined_clone() + MaybeUninit::new(src.vget_unchecked(*source).inlined_clone()); } } } - res + // Safety: res is now fully initialized. + unsafe { res.assume_init() } } /// Creates a new matrix by extracting the given set of columns from `self`. @@ -89,27 +89,30 @@ impl> Matrix { DefaultAllocator: Allocator, { let icols = icols.into_iter(); - let nrows = self.data.shape().0; - let mut res = unsafe { - crate::unimplemented_or_uninitialized_generic!(nrows, Dynamic::new(icols.len())) - }; + let nrows = self.shape_generic().0; + let mut res = Matrix::uninit(nrows, Dynamic::new(icols.len())); for (destination, source) in icols.enumerate() { - res.column_mut(destination).copy_from(&self.column(*source)) + // NOTE: this is basically a copy_frow but wrapping the values insnide of MaybeUninit. + res.column_mut(destination) + .zip_apply(&self.column(*source), |out, e| { + *out = MaybeUninit::new(e.inlined_clone()) + }); } - res + // Safety: res is now fully initialized. + unsafe { res.assume_init() } } } /// # Set rows, columns, and diagonal -impl> Matrix { +impl> Matrix { /// Fills the diagonal of this matrix with the content of the given vector. #[inline] pub fn set_diagonal(&mut self, diag: &Vector) where R: DimMin, - S2: Storage, + S2: RawStorage, ShapeConstraint: DimEq, R2>, { let (nrows, ncols) = self.shape(); @@ -140,7 +143,7 @@ impl> Matrix { #[inline] pub fn set_row(&mut self, i: usize, row: &RowVector) where - S2: Storage, + S2: RawStorage, ShapeConstraint: SameNumberOfColumns, { self.row_mut(i).copy_from(row); @@ -150,7 +153,7 @@ impl> Matrix { #[inline] pub fn set_column(&mut self, i: usize, column: &Vector) where - S2: Storage, + S2: RawStorage, ShapeConstraint: SameNumberOfRows, { self.column_mut(i).copy_from(column); @@ -158,10 +161,21 @@ impl> Matrix { } /// # In-place filling -impl> Matrix { +impl> Matrix { + /// Sets all the elements of this matrix to the value returned by the closure. + #[inline] + pub fn fill_with(&mut self, val: impl Fn() -> T) { + for e in self.iter_mut() { + *e = val() + } + } + /// Sets all the elements of this matrix to `val`. #[inline] - pub fn fill(&mut self, val: T) { + pub fn fill(&mut self, val: T) + where + T: Scalar, + { for e in self.iter_mut() { *e = val.inlined_clone() } @@ -171,7 +185,7 @@ impl> Matrix { #[inline] pub fn fill_with_identity(&mut self) where - T: Zero + One, + T: Scalar + Zero + One, { self.fill(T::zero()); self.fill_diagonal(T::one()); @@ -179,7 +193,10 @@ impl> Matrix { /// Sets all the diagonal elements of this matrix to `val`. #[inline] - pub fn fill_diagonal(&mut self, val: T) { + pub fn fill_diagonal(&mut self, val: T) + where + T: Scalar, + { let (nrows, ncols) = self.shape(); let n = cmp::min(nrows, ncols); @@ -190,7 +207,10 @@ impl> Matrix { /// Sets all the elements of the selected row to `val`. #[inline] - pub fn fill_row(&mut self, i: usize, val: T) { + pub fn fill_row(&mut self, i: usize, val: T) + where + T: Scalar, + { assert!(i < self.nrows(), "Row index out of bounds."); for j in 0..self.ncols() { unsafe { *self.get_unchecked_mut((i, j)) = val.inlined_clone() } @@ -199,7 +219,10 @@ impl> Matrix { /// Sets all the elements of the selected column to `val`. #[inline] - pub fn fill_column(&mut self, j: usize, val: T) { + pub fn fill_column(&mut self, j: usize, val: T) + where + T: Scalar, + { assert!(j < self.ncols(), "Row index out of bounds."); for i in 0..self.nrows() { unsafe { *self.get_unchecked_mut((i, j)) = val.inlined_clone() } @@ -214,7 +237,10 @@ impl> Matrix { /// * If `shift > 1`, then the diagonal and the first `shift - 1` subdiagonals are left /// untouched. #[inline] - pub fn fill_lower_triangle(&mut self, val: T, shift: usize) { + pub fn fill_lower_triangle(&mut self, val: T, shift: usize) + where + T: Scalar, + { for j in 0..self.ncols() { for i in (j + shift)..self.nrows() { unsafe { *self.get_unchecked_mut((i, j)) = val.inlined_clone() } @@ -230,7 +256,10 @@ impl> Matrix { /// * If `shift > 1`, then the diagonal and the first `shift - 1` superdiagonals are left /// untouched. #[inline] - pub fn fill_upper_triangle(&mut self, val: T, shift: usize) { + pub fn fill_upper_triangle(&mut self, val: T, shift: usize) + where + T: Scalar, + { for j in shift..self.ncols() { // TODO: is there a more efficient way to avoid the min ? // (necessary for rectangular matrices) @@ -241,7 +270,7 @@ impl> Matrix { } } -impl> Matrix { +impl> Matrix { /// Copies the upper-triangle of this matrix to its lower-triangular part. /// /// This makes the matrix symmetric. Panics if the matrix is not square. @@ -275,7 +304,7 @@ impl> Matrix { } /// # In-place swapping -impl> Matrix { +impl> Matrix { /// Swaps two rows in-place. #[inline] pub fn swap_rows(&mut self, irow1: usize, irow2: usize) { @@ -335,29 +364,46 @@ impl> Matrix { DefaultAllocator: Reallocator, { let mut m = self.into_owned(); - let (nrows, ncols) = m.data.shape(); + let (nrows, ncols) = m.shape_generic(); let mut offset: usize = 0; let mut target: usize = 0; while offset + target < ncols.value() { if indices.contains(&(target + offset)) { + // Safety: the resulting pointer is within range. + let col_ptr = unsafe { m.data.ptr_mut().add((target + offset) * nrows.value()) }; + // Drop every element in the column we are about to overwrite. + // We use the a similar technique as in `Vec::truncate`. + let s = ptr::slice_from_raw_parts_mut(col_ptr, nrows.value()); + // Safety: we drop the column in-place, which is OK because we will overwrite these + // entries later in the loop, or discard them with the `reallocate_copy` + // afterwards. + unsafe { ptr::drop_in_place(s) }; + offset += 1; } else { unsafe { let ptr_source = m.data.ptr().add((target + offset) * nrows.value()); let ptr_target = m.data.ptr_mut().add(target * nrows.value()); + // Copy the data, overwriting what we dropped. ptr::copy(ptr_source, ptr_target, nrows.value()); target += 1; } } } + // Safety: The new size is smaller than the old size, so + // DefaultAllocator::reallocate_copy will initialize + // every element of the new matrix which can then + // be assumed to be initialized. unsafe { - Matrix::from_data(DefaultAllocator::reallocate_copy( + let new_data = DefaultAllocator::reallocate_copy( nrows, ncols.sub(Dynamic::from_usize(offset)), m.data, - )) + ); + + Matrix::from_data(new_data).assume_init() } } @@ -369,29 +415,44 @@ impl> Matrix { DefaultAllocator: Reallocator, { let mut m = self.into_owned(); - let (nrows, ncols) = m.data.shape(); + let (nrows, ncols) = m.shape_generic(); let mut offset: usize = 0; let mut target: usize = 0; while offset + target < nrows.value() * ncols.value() { if indices.contains(&((target + offset) % nrows.value())) { + // Safety: the resulting pointer is within range. + unsafe { + let elt_ptr = m.data.ptr_mut().add(target + offset); + // Safety: we drop the component in-place, which is OK because we will overwrite these + // entries later in the loop, or discard them with the `reallocate_copy` + // afterwards. + ptr::drop_in_place(elt_ptr) + }; offset += 1; } else { unsafe { let ptr_source = m.data.ptr().add(target + offset); let ptr_target = m.data.ptr_mut().add(target); + // Copy the data, overwriting what we dropped in the previous iterations. ptr::copy(ptr_source, ptr_target, 1); target += 1; } } } + // Safety: The new size is smaller than the old size, so + // DefaultAllocator::reallocate_copy will initialize + // every element of the new matrix which can then + // be assumed to be initialized. unsafe { - Matrix::from_data(DefaultAllocator::reallocate_copy( + let new_data = DefaultAllocator::reallocate_copy( nrows.sub(Dynamic::from_usize(offset / ncols.value())), ncols, m.data, - )) + ); + + Matrix::from_data(new_data).assume_init() } } @@ -432,13 +493,14 @@ impl> Matrix { DefaultAllocator: Reallocator>, { let mut m = self.into_owned(); - let (nrows, ncols) = m.data.shape(); + let (nrows, ncols) = m.shape_generic(); assert!( i + nremove.value() <= ncols.value(), "Column index out of range." ); - if nremove.value() != 0 && i + nremove.value() < ncols.value() { + let need_column_shifts = nremove.value() != 0 && i + nremove.value() < ncols.value(); + if need_column_shifts { // The first `deleted_i * nrows` are left untouched. let copied_value_start = i + nremove.value(); @@ -446,20 +508,35 @@ impl> Matrix { let ptr_in = m.data.ptr().add(copied_value_start * nrows.value()); let ptr_out = m.data.ptr_mut().add(i * nrows.value()); + // Drop all the elements of the columns we are about to overwrite. + // We use the a similar technique as in `Vec::truncate`. + let s = ptr::slice_from_raw_parts_mut(ptr_out, nremove.value() * nrows.value()); + // Safety: we drop the column in-place, which is OK because we will overwrite these + // entries with `ptr::copy` afterward. + ptr::drop_in_place(s); + ptr::copy( ptr_in, ptr_out, (ncols.value() - copied_value_start) * nrows.value(), ); } + } else { + // All the columns to remove are at the end of the buffer. Drop them. + unsafe { + let ptr = m.data.ptr_mut().add(i * nrows.value()); + let s = ptr::slice_from_raw_parts_mut(ptr, nremove.value() * nrows.value()); + ptr::drop_in_place(s) + }; } + // Safety: The new size is smaller than the old size, so + // DefaultAllocator::reallocate_copy will initialize + // every element of the new matrix which can then + // be assumed to be initialized. unsafe { - Matrix::from_data(DefaultAllocator::reallocate_copy( - nrows, - ncols.sub(nremove), - m.data, - )) + let new_data = DefaultAllocator::reallocate_copy(nrows, ncols.sub(nremove), m.data); + Matrix::from_data(new_data).assume_init() } } @@ -511,7 +588,7 @@ impl> Matrix { DefaultAllocator: Reallocator, C>, { let mut m = self.into_owned(); - let (nrows, ncols) = m.data.shape(); + let (nrows, ncols) = m.shape_generic(); assert!( i + nremove.value() <= nrows.value(), "Row index out of range." @@ -520,7 +597,7 @@ impl> Matrix { if nremove.value() != 0 { unsafe { compress_rows( - &mut m.data.as_mut_slice(), + &mut m.as_mut_slice(), nrows.value(), ncols.value(), i, @@ -529,12 +606,13 @@ impl> Matrix { } } + // Safety: The new size is smaller than the old size, so + // DefaultAllocator::reallocate_copy will initialize + // every element of the new matrix which can then + // be assumed to be initialized. unsafe { - Matrix::from_data(DefaultAllocator::reallocate_copy( - nrows.sub(nremove), - ncols, - m.data, - )) + let new_data = DefaultAllocator::reallocate_copy(nrows.sub(nremove), ncols, m.data); + Matrix::from_data(new_data).assume_init() } } } @@ -568,8 +646,13 @@ impl> Matrix { DefaultAllocator: Reallocator>>, { let mut res = unsafe { self.insert_columns_generic_uninitialized(i, Const::) }; - res.fixed_columns_mut::(i).fill(val); - res + res.fixed_columns_mut::(i) + .fill_with(|| MaybeUninit::new(val.inlined_clone())); + + // Safety: the result is now fully initialized. The added columns have + // been initialized by the `fill_with` above, and the rest have + // been initialized by `insert_columns_generic_uninitialized`. + unsafe { res.assume_init() } } /// Inserts `n` columns filled with `val` starting at the `i-th` position. @@ -581,27 +664,33 @@ impl> Matrix { DefaultAllocator: Reallocator, { let mut res = unsafe { self.insert_columns_generic_uninitialized(i, Dynamic::new(n)) }; - res.columns_mut(i, n).fill(val); - res + res.columns_mut(i, n) + .fill_with(|| MaybeUninit::new(val.inlined_clone())); + + // Safety: the result is now fully initialized. The added columns have + // been initialized by the `fill_with` above, and the rest have + // been initialized by `insert_columns_generic_uninitialized`. + unsafe { res.assume_init() } } /// Inserts `ninsert.value()` columns starting at the `i-th` place of this matrix. /// /// # Safety - /// The added column values are not initialized. + /// The output matrix has all its elements initialized except for the the components of the + /// added columns. #[inline] pub unsafe fn insert_columns_generic_uninitialized( self, i: usize, ninsert: D, - ) -> OMatrix> + ) -> UninitMatrix> where D: Dim, C: DimAdd, DefaultAllocator: Reallocator>, { let m = self.into_owned(); - let (nrows, ncols) = m.data.shape(); + let (nrows, ncols) = m.shape_generic(); let mut res = Matrix::from_data(DefaultAllocator::reallocate_copy( nrows, ncols.add(ninsert), @@ -650,8 +739,13 @@ impl> Matrix { DefaultAllocator: Reallocator>, C>, { let mut res = unsafe { self.insert_rows_generic_uninitialized(i, Const::) }; - res.fixed_rows_mut::(i).fill(val); - res + res.fixed_rows_mut::(i) + .fill_with(|| MaybeUninit::new(val.inlined_clone())); + + // Safety: the result is now fully initialized. The added rows have + // been initialized by the `fill_with` above, and the rest have + // been initialized by `insert_rows_generic_uninitialized`. + unsafe { res.assume_init() } } /// Inserts `n` rows filled with `val` starting at the `i-th` position. @@ -663,8 +757,13 @@ impl> Matrix { DefaultAllocator: Reallocator, { let mut res = unsafe { self.insert_rows_generic_uninitialized(i, Dynamic::new(n)) }; - res.rows_mut(i, n).fill(val); - res + res.rows_mut(i, n) + .fill_with(|| MaybeUninit::new(val.inlined_clone())); + + // Safety: the result is now fully initialized. The added rows have + // been initialized by the `fill_with` above, and the rest have + // been initialized by `insert_rows_generic_uninitialized`. + unsafe { res.assume_init() } } /// Inserts `ninsert.value()` rows at the `i-th` place of this matrix. @@ -678,14 +777,14 @@ impl> Matrix { self, i: usize, ninsert: D, - ) -> OMatrix, C> + ) -> UninitMatrix, C> where D: Dim, R: DimAdd, DefaultAllocator: Reallocator, C>, { let m = self.into_owned(); - let (nrows, ncols) = m.data.shape(); + let (nrows, ncols) = m.shape_generic(); let mut res = Matrix::from_data(DefaultAllocator::reallocate_copy( nrows.add(ninsert), ncols, @@ -696,7 +795,7 @@ impl> Matrix { if ninsert.value() != 0 { extend_rows( - &mut res.data.as_mut_slice(), + &mut res.as_mut_slice(), nrows.value(), ncols.value(), i, @@ -731,7 +830,7 @@ impl> Matrix { where DefaultAllocator: Reallocator, { - let ncols = self.data.shape().1; + let ncols = self.shape_generic().1; self.resize_generic(Dynamic::new(new_nrows), ncols, val) } @@ -744,7 +843,7 @@ impl> Matrix { where DefaultAllocator: Reallocator, { - let nrows = self.data.shape().0; + let nrows = self.shape_generic().0; self.resize_generic(nrows, Dynamic::new(new_ncols), val) } @@ -777,16 +876,32 @@ impl> Matrix { DefaultAllocator: Reallocator, { let (nrows, ncols) = self.shape(); - let mut data = self.data.into_owned(); + let mut data = self.into_owned(); if new_nrows.value() == nrows { - let res = unsafe { DefaultAllocator::reallocate_copy(new_nrows, new_ncols, data) }; - let mut res = Matrix::from_data(res); - if new_ncols.value() > ncols { - res.columns_range_mut(ncols..).fill(val); + if new_ncols.value() < ncols { + unsafe { + let num_cols_to_delete = ncols - new_ncols.value(); + let col_ptr = data.data.ptr_mut().add(new_ncols.value() * nrows); + let s = ptr::slice_from_raw_parts_mut(col_ptr, num_cols_to_delete * nrows); + // Safety: drop the elements of the deleted columns. + // these are the elements that will be truncated + // by the `reallocate_copy` afterward. + ptr::drop_in_place(s) + }; } - res + let res = unsafe { DefaultAllocator::reallocate_copy(new_nrows, new_ncols, data.data) }; + let mut res = Matrix::from_data(res); + + if new_ncols.value() > ncols { + res.columns_range_mut(ncols..) + .fill_with(|| MaybeUninit::new(val.inlined_clone())); + } + + // Safety: the result is now fully initialized by `reallocate_copy` and + // `fill_with` (if the output has more columns than the input). + unsafe { res.assume_init() } } else { let mut res; @@ -800,14 +915,14 @@ impl> Matrix { nrows - new_nrows.value(), ); res = Matrix::from_data(DefaultAllocator::reallocate_copy( - new_nrows, new_ncols, data, + new_nrows, new_ncols, data.data, )); } else { res = Matrix::from_data(DefaultAllocator::reallocate_copy( - new_nrows, new_ncols, data, + new_nrows, new_ncols, data.data, )); extend_rows( - &mut res.data.as_mut_slice(), + &mut res.as_mut_slice(), nrows, new_ncols.value(), nrows, @@ -817,15 +932,18 @@ impl> Matrix { } if new_ncols.value() > ncols { - res.columns_range_mut(ncols..).fill(val.inlined_clone()); + res.columns_range_mut(ncols..) + .fill_with(|| MaybeUninit::new(val.inlined_clone())); } if new_nrows.value() > nrows { res.slice_range_mut(nrows.., ..cmp::min(ncols, new_ncols.value())) - .fill(val); + .fill_with(|| MaybeUninit::new(val.inlined_clone())); } - res + // Safety: the result is now fully initialized by `reallocate_copy` and + // `fill_with` (whenever applicable). + unsafe { res.assume_init() } } } @@ -910,12 +1028,8 @@ impl OMatrix { where DefaultAllocator: Reallocator, { - let placeholder = unsafe { - crate::unimplemented_or_uninitialized_generic!(Dynamic::new(0), Dynamic::new(0)) - }; - let old = mem::replace(self, placeholder); - let new = old.resize(new_nrows, new_ncols, val); - let _ = mem::replace(self, new); + // TODO: avoid the clone. + *self = self.clone().resize(new_nrows, new_ncols, val); } } @@ -935,12 +1049,8 @@ where where DefaultAllocator: Reallocator, { - let placeholder = unsafe { - crate::unimplemented_or_uninitialized_generic!(Dynamic::new(0), self.data.shape().1) - }; - let old = mem::replace(self, placeholder); - let new = old.resize_vertically(new_nrows, val); - let _ = mem::replace(self, new); + // TODO: avoid the clone. + *self = self.clone().resize_vertically(new_nrows, val); } } @@ -960,15 +1070,15 @@ where where DefaultAllocator: Reallocator, { - let placeholder = unsafe { - crate::unimplemented_or_uninitialized_generic!(self.data.shape().0, Dynamic::new(0)) - }; - let old = mem::replace(self, placeholder); - let new = old.resize_horizontally(new_ncols, val); - let _ = mem::replace(self, new); + // TODO: avoid the clone. + *self = self.clone().resize_horizontally(new_ncols, val); } } +// Move the elements of `data` in such a way that the matrix with +// the rows `[i, i + nremove[` deleted is represented in a contigous +// way in `data` after this method completes. +// Every deleted element are manually dropped by this method. unsafe fn compress_rows( data: &mut [T], nrows: usize, @@ -978,16 +1088,28 @@ unsafe fn compress_rows( ) { let new_nrows = nrows - nremove; - if new_nrows == 0 || ncols == 0 { - return; // Nothing to do as the output matrix is empty. + if nremove == 0 { + return; // Nothing to remove or drop. } + if new_nrows == 0 || ncols == 0 { + // The output matrix is empty, drop everything. + ptr::drop_in_place(data.as_mut()); + return; + } + + // Safety: because `nremove != 0`, the pointers given to `ptr::copy` + // won’t alias. let ptr_in = data.as_ptr(); let ptr_out = data.as_mut_ptr(); let mut curr_i = i; for k in 0..ncols - 1 { + // Safety: we drop the row elements in-place because we will overwrite these + // entries later with the `ptr::copy`. + let s = ptr::slice_from_raw_parts_mut(ptr_out.add(curr_i), nremove); + ptr::drop_in_place(s); ptr::copy( ptr_in.add(curr_i + (k + 1) * nremove), ptr_out.add(curr_i), @@ -997,7 +1119,13 @@ unsafe fn compress_rows( curr_i += new_nrows; } - // Deal with the last column from which less values have to be copied. + /* + * Deal with the last column from which less values have to be copied. + */ + // Safety: we drop the row elements in-place because we will overwrite these + // entries later with the `ptr::copy`. + let s = ptr::slice_from_raw_parts_mut(ptr_out.add(curr_i), nremove); + ptr::drop_in_place(s); let remaining_len = nrows - i - nremove; ptr::copy( ptr_in.add(nrows * ncols - remaining_len), @@ -1006,15 +1134,9 @@ unsafe fn compress_rows( ); } -// Moves entries of a matrix buffer to make place for `ninsert` emty rows starting at the `i-th` row index. +// Moves entries of a matrix buffer to make place for `ninsert` empty rows starting at the `i-th` row index. // The `data` buffer is assumed to contained at least `(nrows + ninsert) * ncols` elements. -unsafe fn extend_rows( - data: &mut [T], - nrows: usize, - ncols: usize, - i: usize, - ninsert: usize, -) { +unsafe fn extend_rows(data: &mut [T], nrows: usize, ncols: usize, i: usize, ninsert: usize) { let new_nrows = nrows + ninsert; if new_nrows == 0 || ncols == 0 { @@ -1119,7 +1241,7 @@ where R: Dim, S: Extend>, RV: Dim, - SV: Storage, + SV: RawStorage, ShapeConstraint: SameNumberOfRows, { /// Extends the number of columns of a `Matrix` with `Vector`s diff --git a/src/base/indexing.rs b/src/base/indexing.rs index 5107035c..93f41ed3 100644 --- a/src/base/indexing.rs +++ b/src/base/indexing.rs @@ -1,6 +1,6 @@ //! Indexing -use crate::base::storage::{Storage, StorageMut}; +use crate::base::storage::{RawStorage, RawStorageMut}; use crate::base::{ Const, Dim, DimDiff, DimName, DimSub, Dynamic, Matrix, MatrixSlice, MatrixSliceMut, Scalar, U1, }; @@ -310,7 +310,7 @@ fn dimrange_rangetoinclusive_usize() { } /// A helper trait used for indexing operations. -pub trait MatrixIndex<'a, T: Scalar, R: Dim, C: Dim, S: Storage>: Sized { +pub trait MatrixIndex<'a, T, R: Dim, C: Dim, S: RawStorage>: Sized { /// The output type returned by methods. type Output: 'a; @@ -345,7 +345,7 @@ pub trait MatrixIndex<'a, T: Scalar, R: Dim, C: Dim, S: Storage>: Sized } /// A helper trait used for indexing operations. -pub trait MatrixIndexMut<'a, T: Scalar, R: Dim, C: Dim, S: StorageMut>: +pub trait MatrixIndexMut<'a, T, R: Dim, C: Dim, S: RawStorageMut>: MatrixIndex<'a, T, R, C, S> { /// The output type returned by methods. @@ -476,7 +476,7 @@ pub trait MatrixIndexMut<'a, T: Scalar, R: Dim, C: Dim, S: StorageMut>: /// 4, 7, /// 5, 8))); /// ``` -impl> Matrix { +impl> Matrix { /// Produces a view of the data at the given index, or /// `None` if the index is out of bounds. #[inline] @@ -494,7 +494,7 @@ impl> Matrix { #[must_use] pub fn get_mut<'a, I>(&'a mut self, index: I) -> Option where - S: StorageMut, + S: RawStorageMut, I: MatrixIndexMut<'a, T, R, C, S>, { index.get_mut(self) @@ -516,7 +516,7 @@ impl> Matrix { #[inline] pub fn index_mut<'a, I>(&'a mut self, index: I) -> I::OutputMut where - S: StorageMut, + S: RawStorageMut, I: MatrixIndexMut<'a, T, R, C, S>, { index.index_mut(self) @@ -539,7 +539,7 @@ impl> Matrix { #[must_use] pub unsafe fn get_unchecked_mut<'a, I>(&'a mut self, index: I) -> I::OutputMut where - S: StorageMut, + S: RawStorageMut, I: MatrixIndexMut<'a, T, R, C, S>, { index.get_unchecked_mut(self) @@ -553,7 +553,7 @@ where T: Scalar, R: Dim, C: Dim, - S: Storage, + S: RawStorage, { type Output = &'a T; @@ -575,7 +575,7 @@ where T: Scalar, R: Dim, C: Dim, - S: StorageMut, + S: RawStorageMut, { type OutputMut = &'a mut T; @@ -583,7 +583,7 @@ where #[inline(always)] unsafe fn get_unchecked_mut(self, matrix: &'a mut Matrix) -> Self::OutputMut where - S: StorageMut, + S: RawStorageMut, { matrix.data.get_unchecked_linear_mut(self) } @@ -591,12 +591,11 @@ where // EXTRACT A SINGLE ELEMENT BY 2D COORDINATES -impl<'a, T, R, C, S> MatrixIndex<'a, T, R, C, S> for (usize, usize) +impl<'a, T: 'a, R, C, S> MatrixIndex<'a, T, R, C, S> for (usize, usize) where - T: Scalar, R: Dim, C: Dim, - S: Storage, + S: RawStorage, { type Output = &'a T; @@ -604,7 +603,7 @@ where #[inline(always)] fn contained_by(&self, matrix: &Matrix) -> bool { let (rows, cols) = self; - let (nrows, ncols) = matrix.data.shape(); + let (nrows, ncols) = matrix.shape_generic(); DimRange::contained_by(rows, nrows) && DimRange::contained_by(cols, ncols) } @@ -616,12 +615,11 @@ where } } -impl<'a, T, R, C, S> MatrixIndexMut<'a, T, R, C, S> for (usize, usize) +impl<'a, T: 'a, R, C, S> MatrixIndexMut<'a, T, R, C, S> for (usize, usize) where - T: Scalar, R: Dim, C: Dim, - S: StorageMut, + S: RawStorageMut, { type OutputMut = &'a mut T; @@ -629,7 +627,7 @@ where #[inline(always)] unsafe fn get_unchecked_mut(self, matrix: &'a mut Matrix) -> Self::OutputMut where - S: StorageMut, + S: RawStorageMut, { let (row, col) = self; matrix.data.get_unchecked_mut(row, col) @@ -660,7 +658,7 @@ macro_rules! impl_index_pair { T: Scalar, $R: Dim, $C: Dim, - S: Storage, + S: RawStorage, $( $RConstraintType: $RConstraintBound $(<$( $RConstraintBoundParams $( = $REqBound )*),*>)* ,)* $( $CConstraintType: $CConstraintBound $(<$( $CConstraintBoundParams $( = $CEqBound )*),*>)* ),* { @@ -670,7 +668,7 @@ macro_rules! impl_index_pair { #[inline(always)] fn contained_by(&self, matrix: &Matrix) -> bool { let (rows, cols) = self; - let (nrows, ncols) = matrix.data.shape(); + let (nrows, ncols) = matrix.shape_generic(); DimRange::contained_by(rows, nrows) && DimRange::contained_by(cols, ncols) } @@ -680,7 +678,7 @@ macro_rules! impl_index_pair { use crate::base::SliceStorage; let (rows, cols) = self; - let (nrows, ncols) = matrix.data.shape(); + let (nrows, ncols) = matrix.shape_generic(); let data = SliceStorage::new_unchecked(&matrix.data, @@ -696,7 +694,7 @@ macro_rules! impl_index_pair { T: Scalar, $R: Dim, $C: Dim, - S: StorageMut, + S: RawStorageMut, $( $RConstraintType: $RConstraintBound $(<$( $RConstraintBoundParams $( = $REqBound )*),*>)* ,)* $( $CConstraintType: $CConstraintBound $(<$( $CConstraintBoundParams $( = $CEqBound )*),*>)* ),* { @@ -708,7 +706,7 @@ macro_rules! impl_index_pair { use crate::base::SliceStorageMut; let (rows, cols) = self; - let (nrows, ncols) = matrix.data.shape(); + let (nrows, ncols) = matrix.shape_generic(); let data = SliceStorageMut::new_unchecked(&mut matrix.data, diff --git a/src/base/iter.rs b/src/base/iter.rs index 6baeab25..b68e1051 100644 --- a/src/base/iter.rs +++ b/src/base/iter.rs @@ -5,14 +5,14 @@ use std::marker::PhantomData; use std::mem; use crate::base::dimension::{Dim, U1}; -use crate::base::storage::{Storage, StorageMut}; +use crate::base::storage::{RawStorage, RawStorageMut}; use crate::base::{Matrix, MatrixSlice, MatrixSliceMut, Scalar}; macro_rules! iterator { (struct $Name:ident for $Storage:ident.$ptr: ident -> $Ptr:ty, $Ref:ty, $SRef: ty) => { /// An iterator through a dense matrix with arbitrary strides matrix. #[derive(Debug)] - pub struct $Name<'a, T: Scalar, R: Dim, C: Dim, S: 'a + $Storage> { + pub struct $Name<'a, T, R: Dim, C: Dim, S: 'a + $Storage> { ptr: $Ptr, inner_ptr: $Ptr, inner_end: $Ptr, @@ -23,7 +23,7 @@ macro_rules! iterator { // TODO: we need to specialize for the case where the matrix storage is owned (in which // case the iterator is trivial because it does not have any stride). - impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + $Storage> $Name<'a, T, R, C, S> { + impl<'a, T, R: Dim, C: Dim, S: 'a + $Storage> $Name<'a, T, R, C, S> { /// Creates a new iterator for the given matrix storage. pub fn new(storage: $SRef) -> $Name<'a, T, R, C, S> { let shape = storage.shape(); @@ -60,9 +60,7 @@ macro_rules! iterator { } } - impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + $Storage> Iterator - for $Name<'a, T, R, C, S> - { + impl<'a, T, R: Dim, C: Dim, S: 'a + $Storage> Iterator for $Name<'a, T, R, C, S> { type Item = $Ref; #[inline] @@ -117,7 +115,7 @@ macro_rules! iterator { } } - impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + $Storage> DoubleEndedIterator + impl<'a, T, R: Dim, C: Dim, S: 'a + $Storage> DoubleEndedIterator for $Name<'a, T, R, C, S> { #[inline] @@ -157,7 +155,7 @@ macro_rules! iterator { } } - impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + $Storage> ExactSizeIterator + impl<'a, T, R: Dim, C: Dim, S: 'a + $Storage> ExactSizeIterator for $Name<'a, T, R, C, S> { #[inline] @@ -166,15 +164,15 @@ macro_rules! iterator { } } - impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + $Storage> FusedIterator + impl<'a, T, R: Dim, C: Dim, S: 'a + $Storage> FusedIterator for $Name<'a, T, R, C, S> { } }; } -iterator!(struct MatrixIter for Storage.ptr -> *const T, &'a T, &'a S); -iterator!(struct MatrixIterMut for StorageMut.ptr_mut -> *mut T, &'a mut T, &'a mut S); +iterator!(struct MatrixIter for RawStorage.ptr -> *const T, &'a T, &'a S); +iterator!(struct MatrixIterMut for RawStorageMut.ptr_mut -> *mut T, &'a mut T, &'a mut S); /* * @@ -183,18 +181,18 @@ iterator!(struct MatrixIterMut for StorageMut.ptr_mut -> *mut T, &'a mut T, &'a */ #[derive(Clone, Debug)] /// An iterator through the rows of a matrix. -pub struct RowIter<'a, T: Scalar, R: Dim, C: Dim, S: Storage> { +pub struct RowIter<'a, T, R: Dim, C: Dim, S: RawStorage> { mat: &'a Matrix, curr: usize, } -impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + Storage> RowIter<'a, T, R, C, S> { +impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorage> RowIter<'a, T, R, C, S> { pub(crate) fn new(mat: &'a Matrix) -> Self { RowIter { mat, curr: 0 } } } -impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + Storage> Iterator for RowIter<'a, T, R, C, S> { +impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorage> Iterator for RowIter<'a, T, R, C, S> { type Item = MatrixSlice<'a, T, U1, C, S::RStride, S::CStride>; #[inline] @@ -222,7 +220,7 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + Storage> Iterator for RowIt } } -impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + Storage> ExactSizeIterator +impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorage> ExactSizeIterator for RowIter<'a, T, R, C, S> { #[inline] @@ -233,13 +231,13 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + Storage> ExactSizeIterator /// An iterator through the mutable rows of a matrix. #[derive(Debug)] -pub struct RowIterMut<'a, T: Scalar, R: Dim, C: Dim, S: StorageMut> { +pub struct RowIterMut<'a, T, R: Dim, C: Dim, S: RawStorageMut> { mat: *mut Matrix, curr: usize, phantom: PhantomData<&'a mut Matrix>, } -impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + StorageMut> RowIterMut<'a, T, R, C, S> { +impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorageMut> RowIterMut<'a, T, R, C, S> { pub(crate) fn new(mat: &'a mut Matrix) -> Self { RowIterMut { mat, @@ -253,7 +251,7 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + StorageMut> RowIterMut<'a, } } -impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + StorageMut> Iterator +impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorageMut> Iterator for RowIterMut<'a, T, R, C, S> { type Item = MatrixSliceMut<'a, T, U1, C, S::RStride, S::CStride>; @@ -280,7 +278,7 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + StorageMut> Iterator } } -impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + StorageMut> ExactSizeIterator +impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorageMut> ExactSizeIterator for RowIterMut<'a, T, R, C, S> { #[inline] @@ -296,20 +294,18 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + StorageMut> ExactSizeIterat */ #[derive(Clone, Debug)] /// An iterator through the columns of a matrix. -pub struct ColumnIter<'a, T: Scalar, R: Dim, C: Dim, S: Storage> { +pub struct ColumnIter<'a, T, R: Dim, C: Dim, S: RawStorage> { mat: &'a Matrix, curr: usize, } -impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + Storage> ColumnIter<'a, T, R, C, S> { +impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorage> ColumnIter<'a, T, R, C, S> { pub(crate) fn new(mat: &'a Matrix) -> Self { ColumnIter { mat, curr: 0 } } } -impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + Storage> Iterator - for ColumnIter<'a, T, R, C, S> -{ +impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorage> Iterator for ColumnIter<'a, T, R, C, S> { type Item = MatrixSlice<'a, T, R, U1, S::RStride, S::CStride>; #[inline] @@ -337,7 +333,7 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + Storage> Iterator } } -impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + Storage> ExactSizeIterator +impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorage> ExactSizeIterator for ColumnIter<'a, T, R, C, S> { #[inline] @@ -348,13 +344,13 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + Storage> ExactSizeIterator /// An iterator through the mutable columns of a matrix. #[derive(Debug)] -pub struct ColumnIterMut<'a, T: Scalar, R: Dim, C: Dim, S: StorageMut> { +pub struct ColumnIterMut<'a, T, R: Dim, C: Dim, S: RawStorageMut> { mat: *mut Matrix, curr: usize, phantom: PhantomData<&'a mut Matrix>, } -impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + StorageMut> ColumnIterMut<'a, T, R, C, S> { +impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorageMut> ColumnIterMut<'a, T, R, C, S> { pub(crate) fn new(mat: &'a mut Matrix) -> Self { ColumnIterMut { mat, @@ -368,7 +364,7 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + StorageMut> ColumnIterMut<' } } -impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + StorageMut> Iterator +impl<'a, T, R: Dim, C: Dim, S: 'a + RawStorageMut> Iterator for ColumnIterMut<'a, T, R, C, S> { type Item = MatrixSliceMut<'a, T, R, U1, S::RStride, S::CStride>; @@ -395,7 +391,7 @@ impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + StorageMut> Iterator } } -impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + StorageMut> ExactSizeIterator +impl<'a, T: Scalar, R: Dim, C: Dim, S: 'a + RawStorageMut> ExactSizeIterator for ColumnIterMut<'a, T, R, C, S> { #[inline] diff --git a/src/base/matrix.rs b/src/base/matrix.rs index 2f3d9dff..ce5f2f18 100644 --- a/src/base/matrix.rs +++ b/src/base/matrix.rs @@ -25,14 +25,15 @@ use crate::base::dimension::{Dim, DimAdd, DimSum, IsNotStaticOne, U1, U2, U3}; use crate::base::iter::{ ColumnIter, ColumnIterMut, MatrixIter, MatrixIterMut, RowIter, RowIterMut, }; -use crate::base::storage::{ - ContiguousStorage, ContiguousStorageMut, Owned, SameShapeStorage, Storage, StorageMut, -}; +use crate::base::storage::{Owned, RawStorage, RawStorageMut, SameShapeStorage}; use crate::base::{Const, DefaultAllocator, OMatrix, OVector, Scalar, Unit}; -use crate::{ArrayStorage, SMatrix, SimdComplexField}; +use crate::{ArrayStorage, SMatrix, SimdComplexField, Storage, UninitMatrix}; +use crate::storage::IsContiguous; +use crate::uninit::{Init, InitStatus, Uninit}; #[cfg(any(feature = "std", feature = "alloc"))] use crate::{DMatrix, DVector, Dynamic, VecStorage}; +use std::mem::MaybeUninit; /// A square matrix. pub type SquareMatrix = Matrix; @@ -186,7 +187,7 @@ pub struct Matrix { // from_data_statically_unchecked. // Note that it would probably make sense to just have // the type `Matrix`, and have `T, R, C` be associated-types - // of the `Storage` trait. However, because we don't have + // of the `RawStorage` trait. However, because we don't have // specialization, this is not bossible because these `T, R, C` // allows us to desambiguate a lot of configurations. _phantoms: PhantomData<(T, R, C)>, @@ -267,7 +268,7 @@ impl Abomonation for Matrix> matrixcompare_core::Matrix +impl> matrixcompare_core::Matrix for Matrix { fn rows(&self) -> usize { @@ -284,7 +285,7 @@ impl> matrixcompare_core::Matrix< } #[cfg(feature = "compare")] -impl> matrixcompare_core::DenseAccess +impl> matrixcompare_core::DenseAccess for Matrix { fn fetch_single(&self, row: usize, col: usize) -> T { @@ -293,7 +294,7 @@ impl> matrixcompare_core::DenseAc } #[cfg(feature = "bytemuck")] -unsafe impl> bytemuck::Zeroable +unsafe impl> bytemuck::Zeroable for Matrix where S: bytemuck::Zeroable, @@ -301,7 +302,7 @@ where } #[cfg(feature = "bytemuck")] -unsafe impl> bytemuck::Pod for Matrix +unsafe impl> bytemuck::Pod for Matrix where S: bytemuck::Pod, Self: Copy, @@ -410,28 +411,32 @@ impl DVector { } } -impl> Matrix { +impl UninitMatrix +where + DefaultAllocator: Allocator, +{ + /// Assumes a matrix's entries to be initialized. This operation should be near zero-cost. + /// + /// For the similar method that operates on matrix slices, see [`slice_assume_init`]. + /// + /// # Safety + /// The user must make sure that every single entry of the buffer has been initialized, + /// or Undefined Behavior will immediately occur. + #[inline(always)] + pub unsafe fn assume_init(self) -> OMatrix { + OMatrix::from_data(>::assume_init( + self.data, + )) + } +} + +impl> Matrix { /// Creates a new matrix with the given data. #[inline(always)] pub fn from_data(data: S) -> Self { unsafe { Self::from_data_statically_unchecked(data) } } - /// Creates a new uninitialized matrix with the given uninitialized data - pub unsafe fn from_uninitialized_data(data: mem::MaybeUninit) -> mem::MaybeUninit { - let res: Matrix> = Matrix { - data, - _phantoms: PhantomData, - }; - let res: mem::MaybeUninit>> = - mem::MaybeUninit::new(res); - // safety: since we wrap the inner MaybeUninit in an outer MaybeUninit above, the fact that the `data` field is partially-uninitialized is still opaque. - // with s/transmute_copy/transmute/, rustc claims that `MaybeUninit>>` may be of a different size from `MaybeUninit>` - // but MaybeUninit's documentation says "MaybeUninit is guaranteed to have the same size, alignment, and ABI as T", which implies those types should be the same size - let res: mem::MaybeUninit> = mem::transmute_copy(&res); - res - } - /// The shape of this matrix returned as the tuple (number of rows, number of columns). /// /// # Examples: @@ -443,10 +448,17 @@ impl> Matrix { #[inline] #[must_use] pub fn shape(&self) -> (usize, usize) { - let (nrows, ncols) = self.data.shape(); + let (nrows, ncols) = self.shape_generic(); (nrows.value(), ncols.value()) } + /// The shape of this matrix wrapped into their representative types (`Const` or `Dynamic`). + #[inline] + #[must_use] + pub fn shape_generic(&self) -> (R, C) { + self.data.shape() + } + /// The number of rows of this matrix. /// /// # Examples: @@ -573,7 +585,7 @@ impl> Matrix { T: PartialEq, R2: Dim, C2: Dim, - SB: Storage, + SB: RawStorage, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { assert!(self.shape() == other.shape()); @@ -584,6 +596,8 @@ impl> Matrix { #[inline] pub fn into_owned(self) -> OMatrix where + T: Scalar, + S: Storage, DefaultAllocator: Allocator, { Matrix::from_data(self.data.into_owned()) @@ -596,6 +610,8 @@ impl> Matrix { #[inline] pub fn into_owned_sum(self) -> MatrixSum where + T: Scalar, + S: Storage, R2: Dim, C2: Dim, DefaultAllocator: SameShapeAllocator, @@ -621,6 +637,8 @@ impl> Matrix { #[must_use] pub fn clone_owned(&self) -> OMatrix where + T: Scalar, + S: Storage, DefaultAllocator: Allocator, { Matrix::from_data(self.data.clone_owned()) @@ -632,6 +650,8 @@ impl> Matrix { #[must_use] pub fn clone_owned_sum(&self) -> MatrixSum where + T: Scalar, + S: Storage, R2: Dim, C2: Dim, DefaultAllocator: SameShapeAllocator, @@ -641,44 +661,67 @@ impl> Matrix { let nrows: SameShapeR = Dim::from_usize(nrows); let ncols: SameShapeC = Dim::from_usize(ncols); - let mut res: MatrixSum = - unsafe { crate::unimplemented_or_uninitialized_generic!(nrows, ncols) }; + let mut res = Matrix::uninit(nrows, ncols); - // TODO: use copy_from - for j in 0..res.ncols() { - for i in 0..res.nrows() { + unsafe { + // TODO: use copy_from? + for j in 0..res.ncols() { + for i in 0..res.nrows() { + *res.get_unchecked_mut((i, j)) = + MaybeUninit::new(self.get_unchecked((i, j)).inlined_clone()); + } + } + + // SAFETY: the output has been initialized above. + res.assume_init() + } + } + + /// Transposes `self` and store the result into `out`. + #[inline] + fn transpose_to_uninit( + &self, + status: Status, + out: &mut Matrix, + ) where + Status: InitStatus, + T: Scalar, + R2: Dim, + C2: Dim, + SB: RawStorageMut, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, + { + let (nrows, ncols) = self.shape(); + assert!( + (ncols, nrows) == out.shape(), + "Incompatible shape for transposition." + ); + + // TODO: optimize that. + for i in 0..nrows { + for j in 0..ncols { + // Safety: the indices are in range. unsafe { - *res.get_unchecked_mut((i, j)) = self.get_unchecked((i, j)).inlined_clone(); + Status::init( + out.get_unchecked_mut((j, i)), + self.get_unchecked((i, j)).inlined_clone(), + ); } } } - - res } /// Transposes `self` and store the result into `out`. #[inline] pub fn transpose_to(&self, out: &mut Matrix) where + T: Scalar, R2: Dim, C2: Dim, - SB: StorageMut, + SB: RawStorageMut, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { - let (nrows, ncols) = self.shape(); - assert!( - (ncols, nrows) == out.shape(), - "Incompatible shape for transpose-copy." - ); - - // TODO: optimize that. - for i in 0..nrows { - for j in 0..ncols { - unsafe { - *out.get_unchecked_mut((j, i)) = self.get_unchecked((i, j)).inlined_clone(); - } - } - } + self.transpose_to_uninit(Init, out) } /// Transposes `self`. @@ -686,43 +729,43 @@ impl> Matrix { #[must_use = "Did you mean to use transpose_mut()?"] pub fn transpose(&self) -> OMatrix where + T: Scalar, DefaultAllocator: Allocator, { - let (nrows, ncols) = self.data.shape(); + let (nrows, ncols) = self.shape_generic(); - unsafe { - let mut res = crate::unimplemented_or_uninitialized_generic!(ncols, nrows); - self.transpose_to(&mut res); - - res - } + let mut res = Matrix::uninit(ncols, nrows); + self.transpose_to_uninit(Uninit, &mut res); + // Safety: res is now fully initialized. + unsafe { res.assume_init() } } } /// # Elementwise mapping and folding -impl> Matrix { +impl> Matrix { /// Returns a matrix containing the result of `f` applied to each of its entries. #[inline] #[must_use] pub fn map T2>(&self, mut f: F) -> OMatrix where + T: Scalar, DefaultAllocator: Allocator, { - let (nrows, ncols) = self.data.shape(); - - let mut res: OMatrix = - unsafe { crate::unimplemented_or_uninitialized_generic!(nrows, ncols) }; + let (nrows, ncols) = self.shape_generic(); + let mut res = Matrix::uninit(nrows, ncols); for j in 0..ncols.value() { for i in 0..nrows.value() { + // Safety: all indices are in range. unsafe { let a = self.data.get_unchecked(i, j).inlined_clone(); - *res.data.get_unchecked_mut(i, j) = f(a) + *res.data.get_unchecked_mut(i, j) = MaybeUninit::new(f(a)); } } } - res + // Safety: res is now fully initialized. + unsafe { res.assume_init() } } /// Cast the components of `self` to another type. @@ -736,6 +779,7 @@ impl> Matrix { /// ``` pub fn cast(self) -> OMatrix where + T: Scalar, OMatrix: SupersetOf, DefaultAllocator: Allocator, { @@ -755,7 +799,10 @@ impl> Matrix { &self, init_f: impl FnOnce(Option<&T>) -> T2, f: impl FnMut(T2, &T) -> T2, - ) -> T2 { + ) -> T2 + where + T: Scalar, + { let mut it = self.iter(); let init = init_f(it.next()); it.fold(init, f) @@ -770,23 +817,24 @@ impl> Matrix { mut f: F, ) -> OMatrix where + T: Scalar, DefaultAllocator: Allocator, { - let (nrows, ncols) = self.data.shape(); - - let mut res: OMatrix = - unsafe { crate::unimplemented_or_uninitialized_generic!(nrows, ncols) }; + let (nrows, ncols) = self.shape_generic(); + let mut res = Matrix::uninit(nrows, ncols); for j in 0..ncols.value() { for i in 0..nrows.value() { + // Safety: all indices are in range. unsafe { let a = self.data.get_unchecked(i, j).inlined_clone(); - *res.data.get_unchecked_mut(i, j) = f(i, j, a) + *res.data.get_unchecked_mut(i, j) = MaybeUninit::new(f(i, j, a)); } } } - res + // Safety: res is now fully initialized. + unsafe { res.assume_init() } } /// Returns a matrix containing the result of `f` applied to each entries of `self` and @@ -795,16 +843,15 @@ impl> Matrix { #[must_use] pub fn zip_map(&self, rhs: &Matrix, mut f: F) -> OMatrix where + T: Scalar, T2: Scalar, N3: Scalar, - S2: Storage, + S2: RawStorage, F: FnMut(T, T2) -> N3, DefaultAllocator: Allocator, { - let (nrows, ncols) = self.data.shape(); - - let mut res: OMatrix = - unsafe { crate::unimplemented_or_uninitialized_generic!(nrows, ncols) }; + let (nrows, ncols) = self.shape_generic(); + let mut res = Matrix::uninit(nrows, ncols); assert_eq!( (nrows.value(), ncols.value()), @@ -814,15 +861,17 @@ impl> Matrix { for j in 0..ncols.value() { for i in 0..nrows.value() { + // Safety: all indices are in range. unsafe { let a = self.data.get_unchecked(i, j).inlined_clone(); let b = rhs.data.get_unchecked(i, j).inlined_clone(); - *res.data.get_unchecked_mut(i, j) = f(a, b) + *res.data.get_unchecked_mut(i, j) = MaybeUninit::new(f(a, b)) } } } - res + // Safety: res is now fully initialized. + unsafe { res.assume_init() } } /// Returns a matrix containing the result of `f` applied to each entries of `self` and @@ -836,18 +885,17 @@ impl> Matrix { mut f: F, ) -> OMatrix where + T: Scalar, T2: Scalar, N3: Scalar, N4: Scalar, - S2: Storage, - S3: Storage, + S2: RawStorage, + S3: RawStorage, F: FnMut(T, T2, N3) -> N4, DefaultAllocator: Allocator, { - let (nrows, ncols) = self.data.shape(); - - let mut res: OMatrix = - unsafe { crate::unimplemented_or_uninitialized_generic!(nrows, ncols) }; + let (nrows, ncols) = self.shape_generic(); + let mut res = Matrix::uninit(nrows, ncols); assert_eq!( (nrows.value(), ncols.value()), @@ -862,28 +910,34 @@ impl> Matrix { for j in 0..ncols.value() { for i in 0..nrows.value() { + // Safety: all indices are in range. unsafe { let a = self.data.get_unchecked(i, j).inlined_clone(); let b = b.data.get_unchecked(i, j).inlined_clone(); let c = c.data.get_unchecked(i, j).inlined_clone(); - *res.data.get_unchecked_mut(i, j) = f(a, b, c) + *res.data.get_unchecked_mut(i, j) = MaybeUninit::new(f(a, b, c)) } } } - res + // Safety: res is now fully initialized. + unsafe { res.assume_init() } } /// Folds a function `f` on each entry of `self`. #[inline] #[must_use] - pub fn fold(&self, init: Acc, mut f: impl FnMut(Acc, T) -> Acc) -> Acc { - let (nrows, ncols) = self.data.shape(); + pub fn fold(&self, init: Acc, mut f: impl FnMut(Acc, T) -> Acc) -> Acc + where + T: Scalar, + { + let (nrows, ncols) = self.shape_generic(); let mut res = init; for j in 0..ncols.value() { for i in 0..nrows.value() { + // Safety: all indices are in range. unsafe { let a = self.data.get_unchecked(i, j).inlined_clone(); res = f(res, a) @@ -904,13 +958,14 @@ impl> Matrix { mut f: impl FnMut(Acc, T, T2) -> Acc, ) -> Acc where + T: Scalar, T2: Scalar, R2: Dim, C2: Dim, - S2: Storage, + S2: RawStorage, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { - let (nrows, ncols) = self.data.shape(); + let (nrows, ncols) = self.shape_generic(); let mut res = init; @@ -933,11 +988,11 @@ impl> Matrix { res } - /// Replaces each component of `self` by the result of a closure `f` applied on it. + /// Applies a closure `f` to modify each component of `self`. #[inline] - pub fn apply T>(&mut self, mut f: F) + pub fn apply(&mut self, mut f: F) where - S: StorageMut, + S: RawStorageMut, { let (nrows, ncols) = self.shape(); @@ -945,7 +1000,7 @@ impl> Matrix { for i in 0..nrows { unsafe { let e = self.data.get_unchecked_mut(i, j); - *e = f(e.inlined_clone()) + f(e) } } } @@ -957,13 +1012,13 @@ impl> Matrix { pub fn zip_apply( &mut self, rhs: &Matrix, - mut f: impl FnMut(T, T2) -> T, + mut f: impl FnMut(&mut T, T2), ) where - S: StorageMut, + S: RawStorageMut, T2: Scalar, R2: Dim, C2: Dim, - S2: Storage, + S2: RawStorage, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { let (nrows, ncols) = self.shape(); @@ -979,7 +1034,7 @@ impl> Matrix { unsafe { let e = self.data.get_unchecked_mut(i, j); let rhs = rhs.get_unchecked((i, j)).inlined_clone(); - *e = f(e.inlined_clone(), rhs) + f(e, rhs) } } } @@ -992,17 +1047,17 @@ impl> Matrix { &mut self, b: &Matrix, c: &Matrix, - mut f: impl FnMut(T, T2, N3) -> T, + mut f: impl FnMut(&mut T, T2, N3), ) where - S: StorageMut, + S: RawStorageMut, T2: Scalar, R2: Dim, C2: Dim, - S2: Storage, + S2: RawStorage, N3: Scalar, R3: Dim, C3: Dim, - S3: Storage, + S3: RawStorage, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { @@ -1025,7 +1080,7 @@ impl> Matrix { let e = self.data.get_unchecked_mut(i, j); let b = b.get_unchecked((i, j)).inlined_clone(); let c = c.get_unchecked((i, j)).inlined_clone(); - *e = f(e.inlined_clone(), b, c) + f(e, b, c) } } } @@ -1033,7 +1088,7 @@ impl> Matrix { } /// # Iteration on components, rows, and columns -impl> Matrix { +impl> Matrix { /// Iterates through this matrix coordinates in column-major order. /// /// # Examples: @@ -1090,7 +1145,7 @@ impl> Matrix { #[inline] pub fn iter_mut(&mut self) -> MatrixIterMut<'_, T, R, C, S> where - S: StorageMut, + S: RawStorageMut, { MatrixIterMut::new(&mut self.data) } @@ -1113,7 +1168,7 @@ impl> Matrix { #[inline] pub fn row_iter_mut(&mut self) -> RowIterMut<'_, T, R, C, S> where - S: StorageMut, + S: RawStorageMut, { RowIterMut::new(self) } @@ -1136,13 +1191,13 @@ impl> Matrix { #[inline] pub fn column_iter_mut(&mut self) -> ColumnIterMut<'_, T, R, C, S> where - S: StorageMut, + S: RawStorageMut, { ColumnIterMut::new(self) } } -impl> Matrix { +impl> Matrix { /// Returns a mutable pointer to the start of the matrix. /// /// If the matrix is not empty, this pointer is guaranteed to be aligned @@ -1179,7 +1234,10 @@ impl> Matrix { /// /// The components of the slice are assumed to be ordered in column-major order. #[inline] - pub fn copy_from_slice(&mut self, slice: &[T]) { + pub fn copy_from_slice(&mut self, slice: &[T]) + where + T: Scalar, + { let (nrows, ncols) = self.shape(); assert!( @@ -1201,9 +1259,10 @@ impl> Matrix { #[inline] pub fn copy_from(&mut self, other: &Matrix) where + T: Scalar, R2: Dim, C2: Dim, - SB: Storage, + SB: RawStorage, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { assert!( @@ -1224,9 +1283,10 @@ impl> Matrix { #[inline] pub fn tr_copy_from(&mut self, other: &Matrix) where + T: Scalar, R2: Dim, C2: Dim, - SB: Storage, + SB: RawStorage, ShapeConstraint: DimEq + SameNumberOfColumns, { let (nrows, ncols) = self.shape(); @@ -1247,13 +1307,13 @@ impl> Matrix { // TODO: rename `apply` to `apply_mut` and `apply_into` to `apply`? /// Returns `self` with each of its components replaced by the result of a closure `f` applied on it. #[inline] - pub fn apply_into T>(mut self, f: F) -> Self { + pub fn apply_into(mut self, f: F) -> Self { self.apply(f); self } } -impl> Vector { +impl> Vector { /// Gets a reference to the i-th element of this column vector without bound checking. #[inline] #[must_use] @@ -1264,7 +1324,7 @@ impl> Vector { } } -impl> Vector { +impl> Vector { /// Gets a mutable reference to the i-th element of this column vector without bound checking. #[inline] #[must_use] @@ -1275,25 +1335,27 @@ impl> Vector { } } -impl> Matrix { +impl + IsContiguous> Matrix { /// Extracts a slice containing the entire matrix entries ordered column-by-columns. #[inline] #[must_use] pub fn as_slice(&self) -> &[T] { - self.data.as_slice() + // Safety: this is OK thanks to the IsContiguous trait. + unsafe { self.data.as_slice_unchecked() } } } -impl> Matrix { +impl + IsContiguous> Matrix { /// Extracts a mutable slice containing the entire matrix entries ordered column-by-columns. #[inline] #[must_use] pub fn as_mut_slice(&mut self) -> &mut [T] { - self.data.as_mut_slice() + // Safety: this is OK thanks to the IsContiguous trait. + unsafe { self.data.as_mut_slice_unchecked() } } } -impl> Matrix { +impl> Matrix { /// Transposes the square matrix `self` in-place. pub fn transpose_mut(&mut self) { assert!( @@ -1311,14 +1373,18 @@ impl> Matrix { } } -impl> Matrix { +impl> Matrix { /// Takes the adjoint (aka. conjugate-transpose) of `self` and store the result into `out`. #[inline] - pub fn adjoint_to(&self, out: &mut Matrix) - where + fn adjoint_to_uninit( + &self, + status: Status, + out: &mut Matrix, + ) where + Status: InitStatus, R2: Dim, C2: Dim, - SB: StorageMut, + SB: RawStorageMut, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { let (nrows, ncols) = self.shape(); @@ -1330,13 +1396,29 @@ impl> Matrix(&self, out: &mut Matrix) + where + R2: Dim, + C2: Dim, + SB: RawStorageMut, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, + { + self.adjoint_to_uninit(Init, out) + } + /// The adjoint (aka. conjugate-transpose) of `self`. #[inline] #[must_use = "Did you mean to use adjoint_mut()?"] @@ -1344,15 +1426,13 @@ impl> Matrix, { - let (nrows, ncols) = self.data.shape(); + let (nrows, ncols) = self.shape_generic(); - unsafe { - let mut res: OMatrix<_, C, R> = - crate::unimplemented_or_uninitialized_generic!(ncols, nrows); - self.adjoint_to(&mut res); + let mut res = Matrix::uninit(ncols, nrows); + self.adjoint_to_uninit(Uninit, &mut res); - res - } + // Safety: res is now fully initialized. + unsafe { res.assume_init() } } /// Takes the conjugate and transposes `self` and store the result into `out`. @@ -1362,7 +1442,7 @@ impl> Matrix, + SB: RawStorageMut, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { self.adjoint_to(out) @@ -1409,27 +1489,27 @@ impl> Matrix> Matrix { +impl> Matrix { /// The conjugate of the complex matrix `self` computed in-place. #[inline] pub fn conjugate_mut(&mut self) { - self.apply(|e| e.simd_conjugate()) + self.apply(|e| *e = e.simd_conjugate()) } /// Divides each component of the complex matrix `self` by the given real. #[inline] pub fn unscale_mut(&mut self, real: T::SimdRealField) { - self.apply(|e| e.simd_unscale(real)) + self.apply(|e| *e = e.simd_unscale(real)) } /// Multiplies each component of the complex matrix `self` by the given real. #[inline] pub fn scale_mut(&mut self, real: T::SimdRealField) { - self.apply(|e| e.simd_scale(real)) + self.apply(|e| *e = e.simd_scale(real)) } } -impl> Matrix { +impl> Matrix { /// Sets `self` to its adjoint. #[deprecated(note = "Renamed to `self.adjoint_mut()`.")] pub fn conjugate_transform_mut(&mut self) { @@ -1465,7 +1545,7 @@ impl> Matrix { } } -impl> SquareMatrix { +impl> SquareMatrix { /// The diagonal of this matrix. #[inline] #[must_use] @@ -1490,17 +1570,19 @@ impl> SquareMatrix { "Unable to get the diagonal of a non-square matrix." ); - let dim = self.data.shape().0; - let mut res: OVector = - unsafe { crate::unimplemented_or_uninitialized_generic!(dim, Const::<1>) }; + let dim = self.shape_generic().0; + let mut res = Matrix::uninit(dim, Const::<1>); for i in 0..dim.value() { + // Safety: all indices are in range. unsafe { - *res.vget_unchecked_mut(i) = f(self.get_unchecked((i, i)).inlined_clone()); + *res.vget_unchecked_mut(i) = + MaybeUninit::new(f(self.get_unchecked((i, i)).inlined_clone())); } } - res + // Safety: res is now fully initialized. + unsafe { res.assume_init() } } /// Computes a trace of a square matrix, i.e., the sum of its diagonal elements. @@ -1515,7 +1597,7 @@ impl> SquareMatrix { "Cannot compute the trace of non-square matrix." ); - let dim = self.data.shape().0; + let dim = self.shape_generic().0; let mut res = T::zero(); for i in 0..dim.value() { @@ -1563,7 +1645,7 @@ impl> SquareMatrix { } } -impl + IsNotStaticOne, S: Storage> +impl + IsNotStaticOne, S: RawStorage> Matrix { /// Yields the homogeneous matrix for this matrix, i.e., appending an additional dimension and @@ -1580,13 +1662,13 @@ impl + IsNotStaticOne, S: Storage ); let dim = DimSum::::from_usize(self.nrows() + 1); let mut res = OMatrix::identity_generic(dim, dim); - res.generic_slice_mut::((0, 0), self.data.shape()) + res.generic_slice_mut::((0, 0), self.shape_generic()) .copy_from(self); res } } -impl, S: Storage> Vector { +impl, S: RawStorage> Vector { /// Computes the coordinates in projective space of this vector, i.e., appends a `0` to its /// coordinates. #[inline] @@ -1603,7 +1685,7 @@ impl, S: Storage> Vector { #[inline] pub fn from_homogeneous(v: Vector, SB>) -> Option> where - SB: Storage>, + SB: RawStorage>, DefaultAllocator: Allocator, { if v[v.len() - 1].is_zero() { @@ -1615,7 +1697,7 @@ impl, S: Storage> Vector { } } -impl, S: Storage> Vector { +impl, S: RawStorage> Vector { /// Constructs a new vector of higher dimension by appending `element` to the end of `self`. #[inline] #[must_use] @@ -1625,20 +1707,22 @@ impl, S: Storage> Vector { { let len = self.len(); let hnrows = DimSum::::from_usize(len + 1); - let mut res: OVector = - unsafe { crate::unimplemented_or_uninitialized_generic!(hnrows, Const::<1>) }; - res.generic_slice_mut((0, 0), self.data.shape()) - .copy_from(self); - res[(len, 0)] = element; + let mut res = Matrix::uninit(hnrows, Const::<1>); + // This is basically a copy_from except that we warp the copied + // values into MaybeUninit. + res.generic_slice_mut((0, 0), self.shape_generic()) + .zip_apply(self, |out, e| *out = MaybeUninit::new(e)); + res[(len, 0)] = MaybeUninit::new(element); - res + // Safety: res has been fully initialized. + unsafe { res.assume_init() } } } impl AbsDiffEq for Matrix where T: Scalar + AbsDiffEq, - S: Storage, + S: RawStorage, T::Epsilon: Copy, { type Epsilon = T::Epsilon; @@ -1681,7 +1765,7 @@ where impl UlpsEq for Matrix where T: Scalar + UlpsEq, - S: Storage, + S: RawStorage, T::Epsilon: Copy, { #[inline] @@ -1701,7 +1785,7 @@ where impl PartialOrd for Matrix where T: Scalar + PartialOrd, - S: Storage, + S: RawStorage, { #[inline] fn partial_cmp(&self, other: &Self) -> Option { @@ -1793,7 +1877,7 @@ where impl Eq for Matrix where T: Scalar + Eq, - S: Storage, + S: RawStorage, { } @@ -1804,8 +1888,8 @@ where C2: Dim, R: Dim, R2: Dim, - S: Storage, - S2: Storage, + S: RawStorage, + S2: RawStorage, { #[inline] fn eq(&self, right: &Matrix) -> bool { @@ -1818,7 +1902,7 @@ macro_rules! impl_fmt { impl $trait for Matrix where T: Scalar + $trait, - S: Storage, + S: RawStorage, { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { #[cfg(feature = "std")] @@ -1922,7 +2006,7 @@ mod tests { } /// # Cross product -impl> +impl> Matrix { /// The perpendicular product between two 2D column vectors, i.e. `a.x * b.y - a.y * b.x`. @@ -1932,7 +2016,7 @@ impl, + SB: RawStorage, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns + SameNumberOfRows @@ -1962,7 +2046,7 @@ impl, + SB: RawStorage, DefaultAllocator: SameShapeAllocator, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { @@ -1979,8 +2063,7 @@ impl::from_usize(3); let ncols = SameShapeC::::from_usize(1); - let mut res: MatrixCross = - crate::unimplemented_or_uninitialized_generic!(nrows, ncols); + let mut res = Matrix::uninit(nrows, ncols); let ax = self.get_unchecked((0, 0)); let ay = self.get_unchecked((1, 0)); @@ -1990,22 +2073,28 @@ impl::from_usize(1); let ncols = SameShapeC::::from_usize(3); - let mut res: MatrixCross = - crate::unimplemented_or_uninitialized_generic!(nrows, ncols); + let mut res = Matrix::uninit(nrows, ncols); let ax = self.get_unchecked((0, 0)); let ay = self.get_unchecked((0, 1)); @@ -2015,20 +2104,27 @@ impl> Vector { +impl> Vector { /// Computes the matrix `M` such that for all vector `v` we have `M * v == self.cross(&v)`. #[inline] #[must_use] @@ -2073,7 +2169,7 @@ impl> Matrix AbsDiffEq for Unit> where T: Scalar + AbsDiffEq, - S: Storage, + S: RawStorage, T::Epsilon: Copy, { type Epsilon = T::Epsilon; @@ -2115,7 +2211,7 @@ where impl UlpsEq for Unit> where T: Scalar + UlpsEq, - S: Storage, + S: RawStorage, T::Epsilon: Copy, { #[inline] @@ -2134,7 +2230,7 @@ where T: Scalar + Hash, R: Dim, C: Dim, - S: Storage, + S: RawStorage, { fn hash(&self, state: &mut H) { let (nrows, ncols) = self.shape(); diff --git a/src/base/matrix_simba.rs b/src/base/matrix_simba.rs index e0333f45..5c259207 100644 --- a/src/base/matrix_simba.rs +++ b/src/base/matrix_simba.rs @@ -44,7 +44,6 @@ where fn replace(&mut self, i: usize, val: Self::Element) { self.zip_apply(&val, |mut a, b| { a.replace(i, b); - a }) } @@ -52,7 +51,6 @@ where unsafe fn replace_unchecked(&mut self, i: usize, val: Self::Element) { self.zip_apply(&val, |mut a, b| { a.replace_unchecked(i, b); - a }) } diff --git a/src/base/matrix_slice.rs b/src/base/matrix_slice.rs index 29aeb5ec..261d41e2 100644 --- a/src/base/matrix_slice.rs +++ b/src/base/matrix_slice.rs @@ -6,29 +6,29 @@ use crate::base::allocator::Allocator; use crate::base::default_allocator::DefaultAllocator; use crate::base::dimension::{Const, Dim, DimName, Dynamic, IsNotStaticOne, U1}; use crate::base::iter::MatrixIter; -use crate::base::storage::{ContiguousStorage, ContiguousStorageMut, Owned, Storage, StorageMut}; +use crate::base::storage::{IsContiguous, Owned, RawStorage, RawStorageMut, Storage}; use crate::base::{Matrix, Scalar}; macro_rules! slice_storage_impl( ($doc: expr; $Storage: ident as $SRef: ty; $T: ident.$get_addr: ident ($Ptr: ty as $Ref: ty)) => { #[doc = $doc] #[derive(Debug)] - pub struct $T<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> { + pub struct $T<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> { ptr: $Ptr, shape: (R, C), strides: (RStride, CStride), _phantoms: PhantomData<$Ref>, } - unsafe impl<'a, T: Scalar + Send, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send + 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: Scalar + Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync + 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: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> $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 slice without bound checking and from a raw pointer. #[inline] pub unsafe fn from_raw_parts(ptr: $Ptr, @@ -48,7 +48,7 @@ macro_rules! slice_storage_impl( } // Dynamic is arbitrary. It's just to be able to call the constructors with `Slice::` - impl<'a, T: Scalar, R: Dim, C: Dim> $T<'a, T, R, C, Dynamic, Dynamic> { + impl<'a, T, R: Dim, C: Dim> $T<'a, T, R, C, Dynamic, Dynamic> { /// Create a new matrix slice without bound checking. #[inline] pub unsafe fn new_unchecked(storage: $SRef, start: (usize, usize), shape: (R, C)) @@ -78,10 +78,10 @@ macro_rules! slice_storage_impl( } } - impl <'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> + impl <'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> $T<'a, T, R, C, RStride, CStride> where - Self: ContiguousStorage + Self: RawStorage + IsContiguous { /// Extracts the original slice from this storage pub fn into_slice(self) -> &'a [T] { @@ -99,11 +99,11 @@ macro_rules! slice_storage_impl( slice_storage_impl!("A matrix data storage for a matrix slice. Only contains an internal reference \ to another matrix data storage."; - Storage as &'a S; SliceStorage.get_address_unchecked(*const T as &'a T)); + RawStorage as &'a S; SliceStorage.get_address_unchecked(*const T as &'a T)); slice_storage_impl!("A mutable matrix data storage for mutable matrix slice. Only contains an \ internal mutable reference to another matrix data storage."; - StorageMut as &'a mut S; SliceStorageMut.get_address_unchecked_mut(*mut T as &'a mut T) + RawStorageMut as &'a mut S; SliceStorageMut.get_address_unchecked_mut(*mut T as &'a mut T) ); impl<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Copy @@ -128,7 +128,7 @@ impl<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Clone impl<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> SliceStorageMut<'a, T, R, C, RStride, CStride> where - Self: ContiguousStorageMut, + Self: RawStorageMut + IsContiguous, { /// Extracts the original slice from this storage pub fn into_slice_mut(self) -> &'a mut [T] { @@ -144,7 +144,7 @@ where macro_rules! storage_impl( ($($T: ident),* $(,)*) => {$( - unsafe impl<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Storage + unsafe impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> RawStorage for $T<'a, T, R, C, RStride, CStride> { type RStride = RStride; @@ -181,6 +181,21 @@ macro_rules! storage_impl( } } + #[inline] + unsafe fn as_slice_unchecked(&self) -> &[T] { + let (nrows, ncols) = self.shape(); + if nrows.value() != 0 && ncols.value() != 0 { + let sz = self.linear_index(nrows.value() - 1, ncols.value() - 1); + slice::from_raw_parts(self.ptr, sz + 1) + } + else { + slice::from_raw_parts(self.ptr, 0) + } + } + } + + unsafe impl<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Storage + for $T<'a, T, R, C, RStride, CStride> { #[inline] fn into_owned(self) -> Owned where DefaultAllocator: Allocator { @@ -194,25 +209,13 @@ macro_rules! storage_impl( let it = MatrixIter::new(self).cloned(); DefaultAllocator::allocate_from_iterator(nrows, ncols, it) } - - #[inline] - unsafe fn as_slice_unchecked(&self) -> &[T] { - let (nrows, ncols) = self.shape(); - if nrows.value() != 0 && ncols.value() != 0 { - let sz = self.linear_index(nrows.value() - 1, ncols.value() - 1); - slice::from_raw_parts(self.ptr, sz + 1) - } - else { - slice::from_raw_parts(self.ptr, 0) - } - } } )*} ); storage_impl!(SliceStorage, SliceStorageMut); -unsafe impl<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> StorageMut +unsafe impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> RawStorageMut for SliceStorageMut<'a, T, R, C, RStride, CStride> { #[inline] @@ -232,33 +235,22 @@ unsafe impl<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> StorageMu } } -unsafe impl<'a, T: Scalar, R: Dim, CStride: Dim> ContiguousStorage - for SliceStorage<'a, T, R, U1, U1, CStride> -{ -} -unsafe impl<'a, T: Scalar, R: Dim, CStride: Dim> ContiguousStorage - for SliceStorageMut<'a, T, R, U1, U1, CStride> -{ -} -unsafe impl<'a, T: Scalar, R: Dim, CStride: Dim> ContiguousStorageMut +unsafe impl<'a, T, R: Dim, CStride: Dim> IsContiguous for SliceStorage<'a, T, R, U1, U1, CStride> {} +unsafe impl<'a, T, R: Dim, CStride: Dim> IsContiguous for SliceStorageMut<'a, T, R, U1, U1, CStride> { } -unsafe impl<'a, T: Scalar, R: DimName, C: Dim + IsNotStaticOne> ContiguousStorage +unsafe impl<'a, T, R: DimName, C: Dim + IsNotStaticOne> IsContiguous for SliceStorage<'a, T, R, C, U1, R> { } -unsafe impl<'a, T: Scalar, R: DimName, C: Dim + IsNotStaticOne> ContiguousStorage - for SliceStorageMut<'a, T, R, C, U1, R> -{ -} -unsafe impl<'a, T: Scalar, R: DimName, C: Dim + IsNotStaticOne> ContiguousStorageMut +unsafe impl<'a, T, R: DimName, C: Dim + IsNotStaticOne> IsContiguous for SliceStorageMut<'a, T, R, C, U1, R> { } -impl> Matrix { +impl> Matrix { #[inline] fn assert_slice_index( &self, @@ -364,7 +356,7 @@ macro_rules! matrix_slice_impl( pub fn $rows_generic($me: $Me, row_start: usize, nrows: RSlice) -> $MatrixSlice<'_, T, RSlice, C, S::RStride, S::CStride> { - let my_shape = $me.data.shape(); + let my_shape = $me.shape_generic(); $me.assert_slice_index((row_start, 0), (nrows.value(), my_shape.1.value()), (0, 0)); let shape = (nrows, my_shape.1); @@ -382,7 +374,7 @@ macro_rules! matrix_slice_impl( -> $MatrixSlice<'_, T, RSlice, C, Dynamic, S::CStride> where RSlice: Dim { - let my_shape = $me.data.shape(); + let my_shape = $me.shape_generic(); let my_strides = $me.data.strides(); $me.assert_slice_index((row_start, 0), (nrows.value(), my_shape.1.value()), (step, 0)); @@ -452,7 +444,7 @@ macro_rules! matrix_slice_impl( pub fn $columns_generic($me: $Me, first_col: usize, ncols: CSlice) -> $MatrixSlice<'_, T, R, CSlice, S::RStride, S::CStride> { - let my_shape = $me.data.shape(); + let my_shape = $me.shape_generic(); $me.assert_slice_index((0, first_col), (my_shape.0.value(), ncols.value()), (0, 0)); let shape = (my_shape.0, ncols); @@ -469,7 +461,7 @@ macro_rules! matrix_slice_impl( pub fn $columns_generic_with_step($me: $Me, first_col: usize, ncols: CSlice, step: usize) -> $MatrixSlice<'_, T, R, CSlice, S::RStride, Dynamic> { - let my_shape = $me.data.shape(); + let my_shape = $me.shape_generic(); let my_strides = $me.data.strides(); $me.assert_slice_index((0, first_col), (my_shape.0.value(), ncols.value()), (0, step)); @@ -592,7 +584,7 @@ macro_rules! matrix_slice_impl( -> ($MatrixSlice<'_, T, Range1::Size, C, S::RStride, S::CStride>, $MatrixSlice<'_, T, Range2::Size, C, S::RStride, S::CStride>) { - let (nrows, ncols) = $me.data.shape(); + let (nrows, ncols) = $me.shape_generic(); let strides = $me.data.strides(); let start1 = r1.begin(nrows); @@ -628,7 +620,7 @@ macro_rules! matrix_slice_impl( -> ($MatrixSlice<'_, T, R, Range1::Size, S::RStride, S::CStride>, $MatrixSlice<'_, T, R, Range2::Size, S::RStride, S::CStride>) { - let (nrows, ncols) = $me.data.shape(); + let (nrows, ncols) = $me.shape_generic(); let strides = $me.data.strides(); let start1 = r1.begin(ncols); @@ -666,9 +658,9 @@ pub type MatrixSliceMut<'a, T, R, C, RStride = U1, CStride = R> = Matrix>; /// # Slicing based on index and length -impl> Matrix { +impl> Matrix { matrix_slice_impl!( - self: &Self, MatrixSlice, SliceStorage, Storage.get_address_unchecked(), &self.data; + self: &Self, MatrixSlice, SliceStorage, RawStorage.get_address_unchecked(), &self.data; row, row_part, rows, @@ -696,9 +688,9 @@ impl> Matrix { } /// # Mutable slicing based on index and length -impl> Matrix { +impl> Matrix { matrix_slice_impl!( - self: &mut Self, MatrixSliceMut, SliceStorageMut, StorageMut.get_address_unchecked_mut(), &mut self.data; + self: &mut Self, MatrixSliceMut, SliceStorageMut, RawStorageMut.get_address_unchecked_mut(), &mut self.data; row_mut, row_part_mut, rows_mut, @@ -861,7 +853,7 @@ impl SliceRange for RangeInclusive { // TODO: see how much of this overlaps with the general indexing // methods from indexing.rs. -impl> Matrix { +impl> Matrix { /// Slices a sub-matrix containing the rows indexed by the range `rows` and the columns indexed /// by the range `cols`. #[inline] @@ -875,7 +867,7 @@ impl> Matrix { RowRange: SliceRange, ColRange: SliceRange, { - let (nrows, ncols) = self.data.shape(); + let (nrows, ncols) = self.shape_generic(); self.generic_slice( (rows.begin(nrows), cols.begin(ncols)), (rows.size(nrows), cols.size(ncols)), @@ -905,7 +897,7 @@ impl> Matrix { // TODO: see how much of this overlaps with the general indexing // methods from indexing.rs. -impl> Matrix { +impl> Matrix { /// Slices a mutable sub-matrix containing the rows indexed by the range `rows` and the columns /// indexed by the range `cols`. pub fn slice_range_mut( @@ -917,7 +909,7 @@ impl> Matrix { RowRange: SliceRange, ColRange: SliceRange, { - let (nrows, ncols) = self.data.shape(); + let (nrows, ncols) = self.shape_generic(); self.generic_slice_mut( (rows.begin(nrows), cols.begin(ncols)), (rows.size(nrows), cols.size(ncols)), @@ -946,7 +938,6 @@ impl> Matrix { impl<'a, T, R, C, RStride, CStride> From> for MatrixSlice<'a, T, R, C, RStride, CStride> where - T: Scalar, R: Dim, C: Dim, RStride: Dim, diff --git a/src/base/min_max.rs b/src/base/min_max.rs index 83e62d10..3d390194 100644 --- a/src/base/min_max.rs +++ b/src/base/min_max.rs @@ -1,10 +1,10 @@ -use crate::storage::Storage; +use crate::storage::RawStorage; use crate::{ComplexField, Dim, Matrix, Scalar, SimdComplexField, SimdPartialOrd, Vector}; use num::{Signed, Zero}; use simba::simd::SimdSigned; /// # Find the min and max components -impl> Matrix { +impl> Matrix { /// Returns the absolute value of the component with the largest absolute value. /// # Example /// ``` @@ -167,7 +167,7 @@ impl> Matrix { } } -impl> Matrix { +impl> Matrix { /// Computes the index of the matrix component with the largest absolute value. /// /// # Examples: @@ -203,7 +203,7 @@ impl> Matri // TODO: find a way to avoid code duplication just for complex number support. /// # Find the min and max components (vector-specific methods) -impl> Vector { +impl> Vector { /// Computes the index of the vector component with the largest complex or real absolute value. /// /// # Examples: diff --git a/src/base/mod.rs b/src/base/mod.rs index fdfbb5c7..c6279ba3 100644 --- a/src/base/mod.rs +++ b/src/base/mod.rs @@ -33,10 +33,13 @@ mod unit; #[cfg(any(feature = "std", feature = "alloc"))] mod vec_storage; +mod blas_uninit; #[doc(hidden)] pub mod helper; mod interpolation; mod min_max; +/// Mechanisms for working with values that may not be initialized. +pub mod uninit; pub use self::matrix::*; pub use self::norm::*; @@ -50,5 +53,6 @@ pub use self::alias::*; pub use self::alias_slice::*; pub use self::array_storage::*; pub use self::matrix_slice::*; +pub use self::storage::*; #[cfg(any(feature = "std", feature = "alloc"))] pub use self::vec_storage::*; diff --git a/src/base/norm.rs b/src/base/norm.rs index a8548ddd..c138069d 100644 --- a/src/base/norm.rs +++ b/src/base/norm.rs @@ -434,7 +434,7 @@ impl> Matrix { { let n = self.norm(); let le = n.simd_le(min_norm); - self.apply(|e| e.simd_unscale(n).select(le, e)); + self.apply(|e| *e = e.simd_unscale(n).select(le, *e)); SimdOption::new(n, le) } @@ -508,13 +508,8 @@ where /// The i-the canonical basis element. #[inline] fn canonical_basis_element(i: usize) -> Self { - assert!(i < D::dim(), "Index out of bound."); - let mut res = Self::zero(); - unsafe { - *res.data.get_unchecked_linear_mut(i) = T::one(); - } - + res[i] = T::one(); res } diff --git a/src/base/ops.rs b/src/base/ops.rs index f9401cc7..bbeb6d07 100644 --- a/src/base/ops.rs +++ b/src/base/ops.rs @@ -7,20 +7,25 @@ use std::ops::{ use simba::scalar::{ClosedAdd, ClosedDiv, ClosedMul, ClosedNeg, ClosedSub}; use crate::base::allocator::{Allocator, SameShapeAllocator, SameShapeC, SameShapeR}; +use crate::base::blas_uninit::gemm_uninit; use crate::base::constraint::{ AreMultipliable, DimEq, SameNumberOfColumns, SameNumberOfRows, ShapeConstraint, }; use crate::base::dimension::{Dim, DimMul, DimName, DimProd, Dynamic}; -use crate::base::storage::{ContiguousStorageMut, Storage, StorageMut}; +use crate::base::storage::{Storage, StorageMut}; +use crate::base::uninit::Uninit; use crate::base::{DefaultAllocator, Matrix, MatrixSum, OMatrix, Scalar, VectorSlice}; -use crate::SimdComplexField; +use crate::storage::IsContiguous; +use crate::uninit::{Init, InitStatus}; +use crate::{RawStorage, RawStorageMut, SimdComplexField}; +use std::mem::MaybeUninit; /* * * Indexing. * */ -impl> Index for Matrix { +impl> Index for Matrix { type Output = T; #[inline] @@ -30,11 +35,7 @@ impl> Index for Matrix Index<(usize, usize)> for Matrix -where - T: Scalar, - S: Storage, -{ +impl> Index<(usize, usize)> for Matrix { type Output = T; #[inline] @@ -50,7 +51,7 @@ where } // Mutable versions. -impl> IndexMut for Matrix { +impl> IndexMut for Matrix { #[inline] fn index_mut(&mut self, i: usize) -> &mut T { let ij = self.vector_to_matrix_index(i); @@ -58,11 +59,7 @@ impl> IndexMut for Matr } } -impl IndexMut<(usize, usize)> for Matrix -where - T: Scalar, - S: StorageMut, -{ +impl> IndexMut<(usize, usize)> for Matrix { #[inline] fn index_mut(&mut self, ij: (usize, usize)) -> &mut T { let shape = self.shape(); @@ -134,7 +131,7 @@ macro_rules! componentwise_binop_impl( ($Trait: ident, $method: ident, $bound: ident; $TraitAssign: ident, $method_assign: ident, $method_assign_statically_unchecked: ident, $method_assign_statically_unchecked_rhs: ident; - $method_to: ident, $method_to_statically_unchecked: ident) => { + $method_to: ident, $method_to_statically_unchecked_uninit: ident) => { impl> Matrix where T: Scalar + $bound { @@ -147,12 +144,14 @@ macro_rules! componentwise_binop_impl( * */ #[inline] - fn $method_to_statically_unchecked(&self, + fn $method_to_statically_unchecked_uninit(&self, + status: Status, rhs: &Matrix, - out: &mut Matrix) - where SB: Storage, - SC: StorageMut { + out: &mut Matrix) + where Status: InitStatus, + SB: RawStorage, + SC: RawStorageMut { assert_eq!(self.shape(), rhs.shape(), "Matrix addition/subtraction dimensions mismatch."); assert_eq!(self.shape(), out.shape(), "Matrix addition/subtraction output dimensions mismatch."); @@ -164,13 +163,13 @@ macro_rules! componentwise_binop_impl( let arr2 = rhs.data.as_slice_unchecked(); let out = out.data.as_mut_slice_unchecked(); for i in 0 .. arr1.len() { - *out.get_unchecked_mut(i) = arr1.get_unchecked(i).inlined_clone().$method(arr2.get_unchecked(i).inlined_clone()); + Status::init(out.get_unchecked_mut(i), arr1.get_unchecked(i).inlined_clone().$method(arr2.get_unchecked(i).inlined_clone())); } } else { for j in 0 .. self.ncols() { for i in 0 .. self.nrows() { let val = self.get_unchecked((i, j)).inlined_clone().$method(rhs.get_unchecked((i, j)).inlined_clone()); - *out.get_unchecked_mut((i, j)) = val; + Status::init(out.get_unchecked_mut((i, j)), val); } } } @@ -254,7 +253,7 @@ macro_rules! componentwise_binop_impl( SC: StorageMut, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns + SameNumberOfRows + SameNumberOfColumns { - self.$method_to_statically_unchecked(rhs, out) + self.$method_to_statically_unchecked_uninit(Init, rhs, out) } } @@ -320,15 +319,13 @@ macro_rules! componentwise_binop_impl( #[inline] fn $method(self, rhs: &'b Matrix) -> Self::Output { - let mut res = unsafe { - let (nrows, ncols) = self.shape(); - let nrows: SameShapeR = Dim::from_usize(nrows); - let ncols: SameShapeC = Dim::from_usize(ncols); - crate::unimplemented_or_uninitialized_generic!(nrows, ncols) - }; - - self.$method_to_statically_unchecked(rhs, &mut res); - res + let (nrows, ncols) = self.shape(); + let nrows: SameShapeR = Dim::from_usize(nrows); + let ncols: SameShapeC = Dim::from_usize(ncols); + let mut res = Matrix::uninit(nrows, ncols); + self.$method_to_statically_unchecked_uninit(Uninit, rhs, &mut res); + // SAFETY: the output has been initialized above. + unsafe { res.assume_init() } } } @@ -362,10 +359,10 @@ macro_rules! componentwise_binop_impl( componentwise_binop_impl!(Add, add, ClosedAdd; AddAssign, add_assign, add_assign_statically_unchecked, add_assign_statically_unchecked_mut; - add_to, add_to_statically_unchecked); + add_to, add_to_statically_unchecked_uninit); componentwise_binop_impl!(Sub, sub, ClosedSub; SubAssign, sub_assign, sub_assign_statically_unchecked, sub_assign_statically_unchecked_mut; - sub_to, sub_to_statically_unchecked); + sub_to, sub_to_statically_unchecked_uninit); impl iter::Sum for OMatrix where @@ -564,11 +561,12 @@ where #[inline] fn mul(self, rhs: &'b Matrix) -> Self::Output { - let mut res = unsafe { - crate::unimplemented_or_uninitialized_generic!(self.data.shape().0, rhs.data.shape().1) - }; - self.mul_to(rhs, &mut res); - res + let mut res = Matrix::uninit(self.shape_generic().0, rhs.shape_generic().1); + unsafe { + // SAFETY: this is OK because status = Uninit && bevy == 0 + gemm_uninit(Uninit, &mut res, T::one(), self, rhs, T::zero()); + res.assume_init() + } } } @@ -633,7 +631,7 @@ where R2: Dim, T: Scalar + Zero + One + ClosedAdd + ClosedMul, SB: Storage, - SA: ContiguousStorageMut + Clone, + SA: StorageMut + IsContiguous + Clone, // TODO: get rid of the IsContiguous ShapeConstraint: AreMultipliable, DefaultAllocator: Allocator, { @@ -650,7 +648,7 @@ where R2: Dim, T: Scalar + Zero + One + ClosedAdd + ClosedMul, SB: Storage, - SA: ContiguousStorageMut + Clone, + SA: StorageMut + IsContiguous + Clone, // TODO: get rid of the IsContiguous ShapeConstraint: AreMultipliable, // TODO: this is too restrictive. See comments for the non-ref version. DefaultAllocator: Allocator, @@ -676,12 +674,10 @@ where DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { - let mut res = unsafe { - crate::unimplemented_or_uninitialized_generic!(self.data.shape().1, rhs.data.shape().1) - }; - - self.tr_mul_to(rhs, &mut res); - res + let mut res = Matrix::uninit(self.shape_generic().1, rhs.shape_generic().1); + self.xx_mul_to_uninit(Uninit, rhs, &mut res, |a, b| a.dot(b)); + // SAFETY: this is OK because the result is now initialized. + unsafe { res.assume_init() } } /// Equivalent to `self.adjoint() * rhs`. @@ -694,26 +690,26 @@ where DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { - let mut res = unsafe { - crate::unimplemented_or_uninitialized_generic!(self.data.shape().1, rhs.data.shape().1) - }; - - self.ad_mul_to(rhs, &mut res); - res + let mut res = Matrix::uninit(self.shape_generic().1, rhs.shape_generic().1); + self.xx_mul_to_uninit(Uninit, rhs, &mut res, |a, b| a.dotc(b)); + // SAFETY: this is OK because the result is now initialized. + unsafe { res.assume_init() } } #[inline(always)] - fn xx_mul_to( + fn xx_mul_to_uninit( &self, + status: Status, rhs: &Matrix, - out: &mut Matrix, + out: &mut Matrix, dot: impl Fn( &VectorSlice<'_, T, R1, SA::RStride, SA::CStride>, &VectorSlice<'_, T, R2, SB::RStride, SB::CStride>, ) -> T, ) where - SB: Storage, - SC: StorageMut, + Status: InitStatus, + SB: RawStorage, + SC: RawStorageMut, ShapeConstraint: SameNumberOfRows + DimEq + DimEq, { let (nrows1, ncols1) = self.shape(); @@ -742,7 +738,8 @@ where for i in 0..ncols1 { for j in 0..ncols2 { let dot = dot(&self.column(i), &rhs.column(j)); - unsafe { *out.get_unchecked_mut((i, j)) = dot }; + let elt = unsafe { out.get_unchecked_mut((i, j)) }; + Status::init(elt, dot); } } } @@ -759,7 +756,7 @@ where SC: StorageMut, ShapeConstraint: SameNumberOfRows + DimEq + DimEq, { - self.xx_mul_to(rhs, out, |a, b| a.dot(b)) + self.xx_mul_to_uninit(Init, rhs, out, |a, b| a.dot(b)) } /// Equivalent to `self.adjoint() * rhs` but stores the result into `out` to avoid @@ -775,7 +772,7 @@ where SC: StorageMut, ShapeConstraint: SameNumberOfRows + DimEq + DimEq, { - self.xx_mul_to(rhs, out, |a, b| a.dotc(b)) + self.xx_mul_to_uninit(Init, rhs, out, |a, b| a.dotc(b)) } /// Equivalent to `self * rhs` but stores the result into `out` to avoid allocations. @@ -808,34 +805,31 @@ where SB: Storage, DefaultAllocator: Allocator, DimProd>, { - let (nrows1, ncols1) = self.data.shape(); - let (nrows2, ncols2) = rhs.data.shape(); + let (nrows1, ncols1) = self.shape_generic(); + let (nrows2, ncols2) = rhs.shape_generic(); - let mut res = unsafe { - crate::unimplemented_or_uninitialized_generic!(nrows1.mul(nrows2), ncols1.mul(ncols2)) - }; - - { - let mut data_res = res.data.ptr_mut(); + let mut res = Matrix::uninit(nrows1.mul(nrows2), ncols1.mul(ncols2)); + let mut data_res = res.data.ptr_mut(); + unsafe { for j1 in 0..ncols1.value() { for j2 in 0..ncols2.value() { for i1 in 0..nrows1.value() { - unsafe { - let coeff = self.get_unchecked((i1, j1)).inlined_clone(); + let coeff = self.get_unchecked((i1, j1)).inlined_clone(); - for i2 in 0..nrows2.value() { - *data_res = coeff.inlined_clone() - * rhs.get_unchecked((i2, j2)).inlined_clone(); - data_res = data_res.offset(1); - } + for i2 in 0..nrows2.value() { + *data_res = MaybeUninit::new( + coeff.inlined_clone() * rhs.get_unchecked((i2, j2)).inlined_clone(), + ); + data_res = data_res.offset(1); } } } } - } - res + // SAFETY: the result matrix has been initialized by the loop above. + res.assume_init() + } } } diff --git a/src/base/properties.rs b/src/base/properties.rs index 9e250119..091d36ef 100644 --- a/src/base/properties.rs +++ b/src/base/properties.rs @@ -8,8 +8,9 @@ use crate::base::allocator::Allocator; use crate::base::dimension::{Dim, DimMin}; use crate::base::storage::Storage; use crate::base::{DefaultAllocator, Matrix, Scalar, SquareMatrix}; +use crate::RawStorage; -impl> Matrix { +impl> Matrix { /// The total number of elements of this matrix. /// /// # Examples: diff --git a/src/base/scalar.rs b/src/base/scalar.rs index db9e458d..baee6e4f 100644 --- a/src/base/scalar.rs +++ b/src/base/scalar.rs @@ -1,19 +1,10 @@ use std::any::Any; -use std::any::TypeId; use std::fmt::Debug; /// The basic scalar type for all structures of `nalgebra`. /// /// This does not make any assumption on the algebraic properties of `Self`. -pub trait Scalar: Clone + PartialEq + Debug + Any { - #[inline] - /// Tests if `Self` the same as the type `T` - /// - /// Typically used to test of `Self` is a f32 or a f64 with `T::is::()`. - fn is() -> bool { - TypeId::of::() == TypeId::of::() - } - +pub trait Scalar: 'static + Clone + PartialEq + Debug { #[inline(always)] /// Performance hack: Clone doesn't get inlined for Copy types in debug mode, so make it inline anyway. fn inlined_clone(&self) -> Self { diff --git a/src/base/statistics.rs b/src/base/statistics.rs index dbb2231c..ebf694a5 100644 --- a/src/base/statistics.rs +++ b/src/base/statistics.rs @@ -1,11 +1,12 @@ use crate::allocator::Allocator; -use crate::storage::Storage; +use crate::storage::RawStorage; use crate::{Const, DefaultAllocator, Dim, Matrix, OVector, RowOVector, Scalar, VectorSlice, U1}; use num::Zero; use simba::scalar::{ClosedAdd, Field, SupersetOf}; +use std::mem::MaybeUninit; /// # Folding on columns and rows -impl> Matrix { +impl> Matrix { /// Returns a row vector where each element is the result of the application of `f` on the /// corresponding column of the original matrix. #[inline] @@ -17,18 +18,19 @@ impl> Matrix { where DefaultAllocator: Allocator, { - let ncols = self.data.shape().1; - let mut res: RowOVector = - unsafe { crate::unimplemented_or_uninitialized_generic!(Const::<1>, ncols) }; + let ncols = self.shape_generic().1; + let mut res = Matrix::uninit(Const::<1>, ncols); for i in 0..ncols.value() { // TODO: avoid bound checking of column. + // Safety: all indices are in range. unsafe { - *res.get_unchecked_mut((0, i)) = f(self.column(i)); + *res.get_unchecked_mut((0, i)) = MaybeUninit::new(f(self.column(i))); } } - res + // Safety: res is now fully initialized. + unsafe { res.assume_init() } } /// Returns a column vector where each element is the result of the application of `f` on the @@ -44,18 +46,19 @@ impl> Matrix { where DefaultAllocator: Allocator, { - let ncols = self.data.shape().1; - let mut res: OVector = - unsafe { crate::unimplemented_or_uninitialized_generic!(ncols, Const::<1>) }; + let ncols = self.shape_generic().1; + let mut res = Matrix::uninit(ncols, Const::<1>); for i in 0..ncols.value() { // TODO: avoid bound checking of column. + // Safety: all indices are in range. unsafe { - *res.vget_unchecked_mut(i) = f(self.column(i)); + *res.vget_unchecked_mut(i) = MaybeUninit::new(f(self.column(i))); } } - res + // Safety: res is now fully initialized. + unsafe { res.assume_init() } } /// Returns a column vector resulting from the folding of `f` on each column of this matrix. @@ -80,7 +83,7 @@ impl> Matrix { } /// # Common statistics operations -impl> Matrix { +impl> Matrix { /* * * Sum computation. @@ -180,7 +183,7 @@ impl> Matrix { T: ClosedAdd + Zero, DefaultAllocator: Allocator, { - let nrows = self.data.shape().0; + let nrows = self.shape_generic().0; self.compress_columns(OVector::zeros_generic(nrows, Const::<1>), |out, col| { *out += col; }) @@ -283,10 +286,10 @@ impl> Matrix { T: Field + SupersetOf, DefaultAllocator: Allocator, { - let (nrows, ncols) = self.data.shape(); + let (nrows, ncols) = self.shape_generic(); let mut mean = self.column_mean(); - mean.apply(|e| -(e.inlined_clone() * e)); + mean.apply(|e| *e = -(e.inlined_clone() * e.inlined_clone())); let denom = T::one() / crate::convert::<_, T>(ncols.value() as f64); self.compress_columns(mean, |out, col| { @@ -391,7 +394,7 @@ impl> Matrix { T: Field + SupersetOf, DefaultAllocator: Allocator, { - let (nrows, ncols) = self.data.shape(); + let (nrows, ncols) = self.shape_generic(); let denom = T::one() / crate::convert::<_, T>(ncols.value() as f64); self.compress_columns(OVector::zeros_generic(nrows, Const::<1>), |out, col| { out.axpy(denom.inlined_clone(), &col, T::one()) diff --git a/src/base/storage.rs b/src/base/storage.rs index a750904f..76a60ce3 100644 --- a/src/base/storage.rs +++ b/src/base/storage.rs @@ -1,6 +1,5 @@ //! Abstract definition of a matrix data storage. -use std::fmt::Debug; use std::ptr; use crate::base::allocator::{Allocator, SameShapeC, SameShapeR}; @@ -19,24 +18,30 @@ pub type SameShapeStorage = /// The owned data storage that can be allocated from `S`. pub type Owned = >::Buffer; +/// The owned data storage that can be allocated from `S`. +pub type OwnedUninit = >::BufferUninit; + /// The row-stride of the owned data storage for a buffer of dimension `(R, C)`. pub type RStride = - <>::Buffer as Storage>::RStride; + <>::Buffer as RawStorage>::RStride; /// The column-stride of the owned data storage for a buffer of dimension `(R, C)`. pub type CStride = - <>::Buffer as Storage>::CStride; + <>::Buffer as RawStorage>::CStride; /// The trait shared by all matrix data storage. /// /// TODO: doc +/// In generic code, it is recommended use the `Storage` trait bound instead. The `RawStorage` +/// trait bound is generally used by code that needs to work with storages that contains +/// `MaybeUninit` elements. /// /// Note that `Self` must always have a number of elements compatible with the matrix length (given /// by `R` and `C` if they are known at compile-time). For example, implementors of this trait /// should **not** allow the user to modify the size of the underlying buffer with safe methods /// (for example the `VecStorage::data_mut` method is unsafe because the user could change the /// vector's size so that it no longer contains enough elements: this will lead to UB. -pub unsafe trait Storage: Debug + Sized { +pub unsafe trait RawStorage: Sized { /// The static stride of this storage's rows. type RStride: Dim; @@ -121,7 +126,10 @@ pub unsafe trait Storage: Debug + Sized { /// /// Call the safe alternative `matrix.as_slice()` instead. unsafe fn as_slice_unchecked(&self) -> &[T]; +} +/// Trait shared by all matrix data storage that don’t contain any uninitialized elements. +pub unsafe trait Storage: RawStorage { /// Builds a matrix data storage that does not contain any reference. fn into_owned(self) -> Owned where @@ -135,10 +143,14 @@ pub unsafe trait Storage: Debug + Sized { /// Trait implemented by matrix data storage that can provide a mutable access to its elements. /// +/// In generic code, it is recommended use the `StorageMut` trait bound instead. The +/// `RawStorageMut` trait bound is generally used by code that needs to work with storages that +/// contains `MaybeUninit` elements. +/// /// Note that a mutable access does not mean that the matrix owns its data. For example, a mutable /// matrix slice can provide mutable access to its elements even if it does not own its data (it /// contains only an internal reference to them). -pub unsafe trait StorageMut: Storage { +pub unsafe trait RawStorageMut: RawStorage { /// The matrix mutable data pointer. fn ptr_mut(&mut self) -> *mut T; @@ -213,40 +225,29 @@ pub unsafe trait StorageMut: Storage { unsafe fn as_mut_slice_unchecked(&mut self) -> &mut [T]; } -/// A matrix storage that is stored contiguously in memory. -/// -/// The storage requirement means that for any value of `i` in `[0, nrows * ncols - 1]`, the value -/// `.get_unchecked_linear` returns one of the matrix component. This trait is unsafe because -/// failing to comply to this may cause Undefined Behaviors. -pub unsafe trait ContiguousStorage: - Storage +/// Trait shared by all mutable matrix data storage that don’t contain any uninitialized elements. +pub unsafe trait StorageMut: + Storage + RawStorageMut { - /// Converts this data storage to a contiguous slice. - fn as_slice(&self) -> &[T] { - // SAFETY: this is safe because this trait guarantees the fact - // that the data is stored contiguously. - unsafe { self.as_slice_unchecked() } - } } -/// A mutable matrix storage that is stored contiguously in memory. +unsafe impl StorageMut for S +where + R: Dim, + C: Dim, + S: Storage + RawStorageMut, +{ +} + +/// Marker trait indicating that a storage is stored contiguously in memory. /// /// The storage requirement means that for any value of `i` in `[0, nrows * ncols - 1]`, the value /// `.get_unchecked_linear` returns one of the matrix component. This trait is unsafe because /// failing to comply to this may cause Undefined Behaviors. -pub unsafe trait ContiguousStorageMut: - ContiguousStorage + StorageMut -{ - /// Converts this data storage to a contiguous mutable slice. - fn as_mut_slice(&mut self) -> &mut [T] { - // SAFETY: this is safe because this trait guarantees the fact - // that the data is stored contiguously. - unsafe { self.as_mut_slice_unchecked() } - } -} +pub unsafe trait IsContiguous {} /// A matrix storage that can be reshaped in-place. -pub trait ReshapableStorage: Storage +pub trait ReshapableStorage: RawStorage where T: Scalar, R1: Dim, @@ -255,7 +256,7 @@ where C2: Dim, { /// The reshaped storage type. - type Output: Storage; + type Output: RawStorage; /// Reshapes the storage into the output storage type. fn reshape_generic(self, nrows: R2, ncols: C2) -> Self::Output; diff --git a/src/base/swizzle.rs b/src/base/swizzle.rs index 25d6375f..6ed05d81 100644 --- a/src/base/swizzle.rs +++ b/src/base/swizzle.rs @@ -1,5 +1,5 @@ use crate::base::{DimName, Scalar, ToTypenum, Vector, Vector2, Vector3}; -use crate::storage::Storage; +use crate::storage::RawStorage; use typenum::{self, Cmp, Greater}; macro_rules! impl_swizzle { @@ -19,7 +19,7 @@ macro_rules! impl_swizzle { } /// # Swizzling -impl> Vector +impl> Vector where D: DimName + ToTypenum, { diff --git a/src/base/uninit.rs b/src/base/uninit.rs new file mode 100644 index 00000000..92d246df --- /dev/null +++ b/src/base/uninit.rs @@ -0,0 +1,76 @@ +use std::mem::MaybeUninit; + +/// This trait is used to write code that may work on matrices that may or may not +/// be initialized. +/// +/// This trait is used to describe how a value must be accessed to initialize it or +/// to retrieve a reference or mutable reference. Typically, a function accepting +/// both initialized and uninitialized inputs should have a `Status: InitStatus` +/// type parameter. Then the methods of the `Status` can be used to access the element. +/// +/// # Safety +/// This trait must not be implemented outside of this crate. +pub unsafe trait InitStatus: Copy { + /// The type of the values with the initialization status described by `Self`. + type Value; + + /// Initialize the given element. + fn init(out: &mut Self::Value, t: T); + + /// Retrieve a reference to the element, assuming that it is initialized. + /// + /// # Safety + /// This is unsound if the referenced value isn’t initialized. + unsafe fn assume_init_ref(t: &Self::Value) -> &T; + + /// Retrieve a mutable reference to the element, assuming that it is initialized. + /// + /// # Safety + /// This is unsound if the referenced value isn’t initialized. + unsafe fn assume_init_mut(t: &mut Self::Value) -> &mut T; +} + +#[derive(Copy, Clone, Debug, PartialEq, Eq)] +/// A type implementing `InitStatus` indicating that the value is completely initialized. +pub struct Init; +#[derive(Copy, Clone, Debug, PartialEq, Eq)] +/// A type implementing `InitStatus` indicating that the value is completely unitialized. +pub struct Uninit; + +unsafe impl InitStatus for Init { + type Value = T; + + #[inline(always)] + fn init(out: &mut T, t: T) { + *out = t; + } + + #[inline(always)] + unsafe fn assume_init_ref(t: &T) -> &T { + t + } + + #[inline(always)] + unsafe fn assume_init_mut(t: &mut T) -> &mut T { + t + } +} + +unsafe impl InitStatus for Uninit { + type Value = MaybeUninit; + + #[inline(always)] + fn init(out: &mut MaybeUninit, t: T) { + *out = MaybeUninit::new(t); + } + + #[inline(always)] + unsafe fn assume_init_ref(t: &MaybeUninit) -> &T { + std::mem::transmute(t.as_ptr()) // TODO: use t.assume_init_ref() + } + + #[inline(always)] + unsafe fn assume_init_mut(t: &mut MaybeUninit) -> &mut T { + std::mem::transmute(t.as_mut_ptr()) // TODO: use t.assume_init_mut() + } +} diff --git a/src/base/unit.rs b/src/base/unit.rs index fb2f6efe..fa869c09 100644 --- a/src/base/unit.rs +++ b/src/base/unit.rs @@ -10,7 +10,7 @@ use abomonation::Abomonation; use crate::allocator::Allocator; use crate::base::DefaultAllocator; -use crate::storage::Storage; +use crate::storage::RawStorage; use crate::{Dim, Matrix, OMatrix, RealField, Scalar, SimdComplexField, SimdRealField}; /// A wrapper that ensures the underlying algebraic entity has a unit norm. @@ -116,7 +116,7 @@ where T: Scalar + PartialEq, R: Dim, C: Dim, - S: Storage, + S: RawStorage, { #[inline] fn eq(&self, rhs: &Self) -> bool { @@ -129,7 +129,7 @@ where T: Scalar + Eq, R: Dim, C: Dim, - S: Storage, + S: RawStorage, { } diff --git a/src/base/vec_storage.rs b/src/base/vec_storage.rs index 20cc5171..bf73661d 100644 --- a/src/base/vec_storage.rs +++ b/src/base/vec_storage.rs @@ -8,9 +8,7 @@ use crate::base::allocator::Allocator; use crate::base::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::base::default_allocator::DefaultAllocator; use crate::base::dimension::{Dim, DimName, Dynamic, U1}; -use crate::base::storage::{ - ContiguousStorage, ContiguousStorageMut, Owned, ReshapableStorage, Storage, StorageMut, -}; +use crate::base::storage::{IsContiguous, Owned, RawStorage, RawStorageMut, ReshapableStorage}; use crate::base::{Scalar, Vector}; #[cfg(feature = "serde-serialize-no-std")] @@ -19,12 +17,14 @@ use serde::{ ser::{Serialize, Serializer}, }; +use crate::Storage; #[cfg(feature = "abomonation-serialize")] use abomonation::Abomonation; +use std::mem::MaybeUninit; /* * - * Storage. + * RawStorage. * */ /// A Vec-based matrix data storage. It may be dynamically-sized. @@ -113,21 +113,49 @@ impl VecStorage { /// Resizes the underlying mutable data storage and unwraps it. /// /// # Safety - /// If `sz` is larger than the current size, additional elements are uninitialized. - /// If `sz` is smaller than the current size, additional elements are truncated. + /// - If `sz` is larger than the current size, additional elements are uninitialized. + /// - If `sz` is smaller than the current size, additional elements are truncated but **not** dropped. + /// It is the responsibility of the caller of this method to drop these elements. #[inline] - pub unsafe fn resize(mut self, sz: usize) -> Vec { + pub unsafe fn resize(mut self, sz: usize) -> Vec> { let len = self.len(); - if sz < len { + let new_data = if sz < len { + // Use `set_len` instead of `truncate` because we don’t want to + // drop the removed elements (it’s the caller’s responsibility). self.data.set_len(sz); self.data.shrink_to_fit(); + + // Safety: + // - MaybeUninit has the same alignment and layout as T. + // - The length and capacity come from a valid vector. + Vec::from_raw_parts( + self.data.as_mut_ptr() as *mut MaybeUninit, + self.data.len(), + self.data.capacity(), + ) } else { self.data.reserve_exact(sz - len); - self.data.set_len(sz); - } - self.data + // Safety: + // - MaybeUninit has the same alignment and layout as T. + // - The length and capacity come from a valid vector. + let mut new_data = Vec::from_raw_parts( + self.data.as_mut_ptr() as *mut MaybeUninit, + self.data.len(), + self.data.capacity(), + ); + + // Safety: we can set the length here because MaybeUninit is always assumed + // to be initialized. + new_data.set_len(sz); + new_data + }; + + // Avoid double-free by forgetting `self` because its data buffer has + // been transfered to `new_data`. + std::mem::forget(self); + new_data } /// The number of elements on the underlying vector. @@ -143,6 +171,18 @@ impl VecStorage { pub fn is_empty(&self) -> bool { self.len() == 0 } + + /// A slice containing all the components stored in this storage in column-major order. + #[inline] + pub fn as_slice(&self) -> &[T] { + &self.data[..] + } + + /// A mutable slice containing all the components stored in this storage in column-major order. + #[inline] + pub fn as_mut_slice(&mut self) -> &mut [T] { + &mut self.data[..] + } } impl From> for Vec { @@ -157,10 +197,7 @@ impl From> for Vec { * Dynamic − Dynamic * */ -unsafe impl Storage for VecStorage -where - DefaultAllocator: Allocator, -{ +unsafe impl RawStorage for VecStorage { type RStride = U1; type CStride = Dynamic; @@ -184,6 +221,16 @@ where true } + #[inline] + unsafe fn as_slice_unchecked(&self) -> &[T] { + &self.data + } +} + +unsafe impl Storage for VecStorage +where + DefaultAllocator: Allocator, +{ #[inline] fn into_owned(self) -> Owned where @@ -199,17 +246,9 @@ where { self.clone() } - - #[inline] - unsafe fn as_slice_unchecked(&self) -> &[T] { - &self.data - } } -unsafe impl Storage for VecStorage -where - DefaultAllocator: Allocator, -{ +unsafe impl RawStorage for VecStorage { type RStride = U1; type CStride = R; @@ -233,6 +272,16 @@ where true } + #[inline] + unsafe fn as_slice_unchecked(&self) -> &[T] { + &self.data + } +} + +unsafe impl Storage for VecStorage +where + DefaultAllocator: Allocator, +{ #[inline] fn into_owned(self) -> Owned where @@ -248,22 +297,14 @@ where { self.clone() } - - #[inline] - unsafe fn as_slice_unchecked(&self) -> &[T] { - &self.data - } } /* * - * StorageMut, ContiguousStorage. + * RawStorageMut, ContiguousStorage. * */ -unsafe impl StorageMut for VecStorage -where - DefaultAllocator: Allocator, -{ +unsafe impl RawStorageMut for VecStorage { #[inline] fn ptr_mut(&mut self) -> *mut T { self.data.as_mut_ptr() @@ -275,15 +316,7 @@ where } } -unsafe impl ContiguousStorage for VecStorage where - DefaultAllocator: Allocator -{ -} - -unsafe impl ContiguousStorageMut for VecStorage where - DefaultAllocator: Allocator -{ -} +unsafe impl IsContiguous for VecStorage {} impl ReshapableStorage for VecStorage where @@ -321,10 +354,7 @@ where } } -unsafe impl StorageMut for VecStorage -where - DefaultAllocator: Allocator, -{ +unsafe impl RawStorageMut for VecStorage { #[inline] fn ptr_mut(&mut self) -> *mut T { self.data.as_mut_ptr() @@ -387,16 +417,6 @@ impl Abomonation for VecStorage { } } -unsafe impl ContiguousStorage for VecStorage where - DefaultAllocator: Allocator -{ -} - -unsafe impl ContiguousStorageMut for VecStorage where - DefaultAllocator: Allocator -{ -} - impl Extend for VecStorage { /// Extends the number of columns of the `VecStorage` with elements /// from the given iterator. @@ -431,7 +451,7 @@ where T: Scalar, R: Dim, RV: Dim, - SV: Storage, + SV: RawStorage, ShapeConstraint: SameNumberOfRows, { /// Extends the number of columns of the `VecStorage` with vectors diff --git a/src/geometry/point.rs b/src/geometry/point.rs index 49138028..098b5c2a 100644 --- a/src/geometry/point.rs +++ b/src/geometry/point.rs @@ -18,6 +18,7 @@ use crate::base::allocator::Allocator; use crate::base::dimension::{DimName, DimNameAdd, DimNameSum, U1}; use crate::base::iter::{MatrixIter, MatrixIterMut}; use crate::base::{Const, DefaultAllocator, OVector, Scalar}; +use std::mem::MaybeUninit; /// A point in an euclidean space. /// @@ -162,16 +163,16 @@ where /// ``` /// # use nalgebra::{Point2, Point3}; /// let mut p = Point2::new(1.0, 2.0); - /// p.apply(|e| e * 10.0); + /// p.apply(|e| *e = *e * 10.0); /// assert_eq!(p, Point2::new(10.0, 20.0)); /// /// // This works in any dimension. /// let mut p = Point3::new(1.0, 2.0, 3.0); - /// p.apply(|e| e * 10.0); + /// p.apply(|e| *e = *e * 10.0); /// assert_eq!(p, Point3::new(10.0, 20.0, 30.0)); /// ``` #[inline] - pub fn apply T>(&mut self, f: F) { + pub fn apply(&mut self, f: F) { self.coords.apply(f) } @@ -198,17 +199,20 @@ where D: DimNameAdd, DefaultAllocator: Allocator>, { - let mut res = unsafe { - crate::unimplemented_or_uninitialized_generic!( - as DimName>::name(), - Const::<1> - ) - }; - res.generic_slice_mut((0, 0), (D::name(), Const::<1>)) - .copy_from(&self.coords); - res[(D::dim(), 0)] = T::one(); + // TODO: this is mostly a copy-past from Vector::push. + // But we can’t use Vector::push because of the DimAdd bound + // (which we don’t use because we use DimNameAdd). + // We should find a way to re-use Vector::push. + let len = self.len(); + let mut res = crate::Matrix::uninit(DimNameSum::::name(), Const::<1>); + // This is basically a copy_from except that we warp the copied + // values into MaybeUninit. + res.generic_slice_mut((0, 0), self.coords.shape_generic()) + .zip_apply(&self.coords, |out, e| *out = MaybeUninit::new(e)); + res[(len, 0)] = MaybeUninit::new(T::one()); - res + // Safety: res has been fully initialized. + unsafe { res.assume_init() } } /// Creates a new point with the given coordinates. diff --git a/src/geometry/point_construction.rs b/src/geometry/point_construction.rs index 0ffbf4d8..d2393146 100644 --- a/src/geometry/point_construction.rs +++ b/src/geometry/point_construction.rs @@ -24,15 +24,6 @@ impl OPoint where DefaultAllocator: Allocator, { - /// Creates a new point with uninitialized coordinates. - #[inline] - pub unsafe fn new_uninitialized() -> Self { - Self::from(crate::unimplemented_or_uninitialized_generic!( - D::name(), - Const::<1> - )) - } - /// Creates a new point with all coordinates equal to zero. /// /// # Example diff --git a/src/lib.rs b/src/lib.rs index e21f0709..5fc38070 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -73,11 +73,11 @@ an optimized set of tools for computer graphics and physics. Those features incl #![allow(unused_variables, unused_mut)] #![deny( + missing_docs, nonstandard_style, unused_parens, unused_qualifications, unused_results, - missing_docs, rust_2018_idioms, rust_2018_compatibility, future_incompatible, @@ -88,7 +88,6 @@ an optimized set of tools for computer graphics and physics. Those features incl html_root_url = "https://docs.rs/nalgebra/0.25.0" )] #![cfg_attr(not(feature = "std"), no_std)] -#![cfg_attr(feature = "no_unsound_assume_init", allow(unreachable_code))] #[cfg(feature = "rand-no-std")] extern crate rand_package as rand; diff --git a/src/linalg/balancing.rs b/src/linalg/balancing.rs index f4f8b659..15679e2b 100644 --- a/src/linalg/balancing.rs +++ b/src/linalg/balancing.rs @@ -5,7 +5,6 @@ use std::ops::{DivAssign, MulAssign}; use crate::allocator::Allocator; use crate::base::dimension::Dim; -use crate::base::storage::Storage; use crate::base::{Const, DefaultAllocator, OMatrix, OVector}; /// Applies in-place a modified Parlett and Reinsch matrix balancing with 2-norm to the matrix and returns @@ -18,7 +17,7 @@ where { assert!(matrix.is_square(), "Unable to balance a non-square matrix."); - let dim = matrix.data.shape().0; + let dim = matrix.shape_generic().0; let radix: T = crate::convert(2.0f64); let mut d = OVector::from_element_generic(dim, Const::<1>, T::one()); diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs index 6a462988..e269b4a0 100644 --- a/src/linalg/bidiagonal.rs +++ b/src/linalg/bidiagonal.rs @@ -4,11 +4,11 @@ use serde::{Deserialize, Serialize}; use crate::allocator::Allocator; use crate::base::{DefaultAllocator, Matrix, OMatrix, OVector, Unit}; use crate::dimension::{Const, Dim, DimDiff, DimMin, DimMinimum, DimSub, U1}; -use crate::storage::Storage; use simba::scalar::ComplexField; use crate::geometry::Reflection; use crate::linalg::householder; +use std::mem::MaybeUninit; /// The bidiagonalization of a general matrix. #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] @@ -73,7 +73,7 @@ where { /// Computes the Bidiagonal decomposition using householder reflections. pub fn new(mut matrix: OMatrix) -> Self { - let (nrows, ncols) = matrix.data.shape(); + let (nrows, ncols) = matrix.shape_generic(); let min_nrows_ncols = nrows.min(ncols); let dim = min_nrows_ncols.value(); assert!( @@ -81,68 +81,65 @@ where "Cannot compute the bidiagonalization of an empty matrix." ); - let mut diagonal = - unsafe { crate::unimplemented_or_uninitialized_generic!(min_nrows_ncols, Const::<1>) }; - let mut off_diagonal = unsafe { - crate::unimplemented_or_uninitialized_generic!( - min_nrows_ncols.sub(Const::<1>), - Const::<1> - ) - }; - let mut axis_packed = - unsafe { crate::unimplemented_or_uninitialized_generic!(ncols, Const::<1>) }; - let mut work = unsafe { crate::unimplemented_or_uninitialized_generic!(nrows, Const::<1>) }; + let mut diagonal = Matrix::uninit(min_nrows_ncols, Const::<1>); + let mut off_diagonal = Matrix::uninit(min_nrows_ncols.sub(Const::<1>), Const::<1>); + let mut axis_packed = Matrix::zeros_generic(ncols, Const::<1>); + let mut work = Matrix::zeros_generic(nrows, Const::<1>); let upper_diagonal = nrows.value() >= ncols.value(); if upper_diagonal { for ite in 0..dim - 1 { - householder::clear_column_unchecked(&mut matrix, &mut diagonal[ite], ite, 0, None); - householder::clear_row_unchecked( + diagonal[ite] = MaybeUninit::new(householder::clear_column_unchecked( + &mut matrix, + ite, + 0, + None, + )); + off_diagonal[ite] = MaybeUninit::new(householder::clear_row_unchecked( &mut matrix, - &mut off_diagonal[ite], &mut axis_packed, &mut work, ite, 1, - ); + )); } - householder::clear_column_unchecked( + diagonal[dim - 1] = MaybeUninit::new(householder::clear_column_unchecked( &mut matrix, - &mut diagonal[dim - 1], dim - 1, 0, None, - ); + )); } else { for ite in 0..dim - 1 { - householder::clear_row_unchecked( + diagonal[ite] = MaybeUninit::new(householder::clear_row_unchecked( &mut matrix, - &mut diagonal[ite], &mut axis_packed, &mut work, ite, 0, - ); - householder::clear_column_unchecked( + )); + off_diagonal[ite] = MaybeUninit::new(householder::clear_column_unchecked( &mut matrix, - &mut off_diagonal[ite], ite, 1, None, - ); + )); } - householder::clear_row_unchecked( + diagonal[dim - 1] = MaybeUninit::new(householder::clear_row_unchecked( &mut matrix, - &mut diagonal[dim - 1], &mut axis_packed, &mut work, dim - 1, 0, - ); + )); } + // Safety: diagonal and off_diagonal have been fully initialized. + let (diagonal, off_diagonal) = + unsafe { (diagonal.assume_init(), off_diagonal.assume_init()) }; + Bidiagonal { uv: matrix, diagonal, @@ -194,7 +191,7 @@ where where DefaultAllocator: Allocator, DimMinimum>, { - let (nrows, ncols) = self.uv.data.shape(); + let (nrows, ncols) = self.uv.shape_generic(); let d = nrows.min(ncols); let mut res = OMatrix::identity_generic(d, d); @@ -214,7 +211,7 @@ where where DefaultAllocator: Allocator>, { - let (nrows, ncols) = self.uv.data.shape(); + let (nrows, ncols) = self.uv.shape_generic(); let mut res = Matrix::identity_generic(nrows, nrows.min(ncols)); let dim = self.diagonal.len(); @@ -245,14 +242,12 @@ where where DefaultAllocator: Allocator, C>, { - let (nrows, ncols) = self.uv.data.shape(); + let (nrows, ncols) = self.uv.shape_generic(); let min_nrows_ncols = nrows.min(ncols); let mut res = Matrix::identity_generic(min_nrows_ncols, ncols); - let mut work = - unsafe { crate::unimplemented_or_uninitialized_generic!(min_nrows_ncols, Const::<1>) }; - let mut axis_packed = - unsafe { crate::unimplemented_or_uninitialized_generic!(ncols, Const::<1>) }; + let mut work = Matrix::zeros_generic(min_nrows_ncols, Const::<1>); + let mut axis_packed = Matrix::zeros_generic(ncols, Const::<1>); let shift = self.axis_shift().1; @@ -353,7 +348,7 @@ where // assert!(self.uv.is_square(), "Bidiagonal inverse: unable to compute the inverse of a non-square matrix."); // // // TODO: is there a less naive method ? -// let (nrows, ncols) = self.uv.data.shape(); +// let (nrows, ncols) = self.uv.shape_generic(); // let mut res = OMatrix::identity_generic(nrows, ncols); // self.solve_mut(&mut res); // res diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index f66fb42f..47939311 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -136,7 +136,7 @@ where /// Computes the inverse of the decomposed matrix. #[must_use] pub fn inverse(&self) -> OMatrix { - let shape = self.chol.data.shape(); + let shape = self.chol.shape_generic(); let mut res = OMatrix::identity_generic(shape.0, shape.1); self.solve_mut(&mut res); @@ -237,12 +237,11 @@ where assert!(j < n, "j needs to be within the bound of the new matrix."); // loads the data into a new matrix with an additional jth row/column - let mut chol = unsafe { - crate::unimplemented_or_uninitialized_generic!( - self.chol.data.shape().0.add(Const::<1>), - self.chol.data.shape().1.add(Const::<1>) - ) - }; + // TODO: would it be worth it to avoid the zero-initialization? + let mut chol = Matrix::zeros_generic( + self.chol.shape_generic().0.add(Const::<1>), + self.chol.shape_generic().1.add(Const::<1>), + ); chol.slice_range_mut(..j, ..j) .copy_from(&self.chol.slice_range(..j, ..j)); chol.slice_range_mut(..j, j + 1..) @@ -303,12 +302,11 @@ where assert!(j < n, "j needs to be within the bound of the matrix."); // loads the data into a new matrix except for the jth row/column - let mut chol = unsafe { - crate::unimplemented_or_uninitialized_generic!( - self.chol.data.shape().0.sub(Const::<1>), - self.chol.data.shape().1.sub(Const::<1>) - ) - }; + // TODO: would it be worth it to avoid this zero initialization? + let mut chol = Matrix::zeros_generic( + self.chol.shape_generic().0.sub(Const::<1>), + self.chol.shape_generic().1.sub(Const::<1>), + ); chol.slice_range_mut(..j, ..j) .copy_from(&self.chol.slice_range(..j, ..j)); chol.slice_range_mut(..j, j..) diff --git a/src/linalg/col_piv_qr.rs b/src/linalg/col_piv_qr.rs index fcd0f376..f5c61336 100644 --- a/src/linalg/col_piv_qr.rs +++ b/src/linalg/col_piv_qr.rs @@ -6,11 +6,12 @@ use crate::allocator::{Allocator, Reallocator}; use crate::base::{Const, DefaultAllocator, Matrix, OMatrix, OVector, Unit}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::dimension::{Dim, DimMin, DimMinimum}; -use crate::storage::{Storage, StorageMut}; +use crate::storage::StorageMut; use crate::ComplexField; use crate::geometry::Reflection; use crate::linalg::{householder, PermutationSequence}; +use std::mem::MaybeUninit; /// The QR decomposition (with column pivoting) of a general matrix. #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] @@ -62,30 +63,33 @@ where { /// Computes the `ColPivQR` decomposition using householder reflections. pub fn new(mut matrix: OMatrix) -> Self { - let (nrows, ncols) = matrix.data.shape(); + let (nrows, ncols) = matrix.shape_generic(); let min_nrows_ncols = nrows.min(ncols); let mut p = PermutationSequence::identity_generic(min_nrows_ncols); - let mut diag = - unsafe { crate::unimplemented_or_uninitialized_generic!(min_nrows_ncols, Const::<1>) }; - if min_nrows_ncols.value() == 0 { return ColPivQR { col_piv_qr: matrix, p, - diag, + diag: Matrix::zeros_generic(min_nrows_ncols, Const::<1>), }; } + let mut diag = Matrix::uninit(min_nrows_ncols, Const::<1>); + for i in 0..min_nrows_ncols.value() { let piv = matrix.slice_range(i.., i..).icamax_full(); let col_piv = piv.1 + i; matrix.swap_columns(i, col_piv); p.append_permutation(i, col_piv); - householder::clear_column_unchecked(&mut matrix, &mut diag[i], i, 0, None); + diag[i] = + MaybeUninit::new(householder::clear_column_unchecked(&mut matrix, i, 0, None)); } + // Safety: diag is now fully initialized. + let diag = unsafe { diag.assume_init() }; + ColPivQR { col_piv_qr: matrix, p, @@ -100,7 +104,7 @@ where where DefaultAllocator: Allocator, C>, { - let (nrows, ncols) = self.col_piv_qr.data.shape(); + let (nrows, ncols) = self.col_piv_qr.shape_generic(); let mut res = self .col_piv_qr .rows_generic(0, nrows.min(ncols)) @@ -117,7 +121,7 @@ where where DefaultAllocator: Reallocator, C>, { - let (nrows, ncols) = self.col_piv_qr.data.shape(); + let (nrows, ncols) = self.col_piv_qr.shape_generic(); let mut res = self .col_piv_qr .resize_generic(nrows.min(ncols), ncols, T::zero()); @@ -132,7 +136,7 @@ where where DefaultAllocator: Allocator>, { - let (nrows, ncols) = self.col_piv_qr.data.shape(); + let (nrows, ncols) = self.col_piv_qr.shape_generic(); // NOTE: we could build the identity matrix and call q_mul on it. // Instead we don't so that we take in account the matrix sparseness. @@ -295,7 +299,7 @@ where ); // TODO: is there a less naive method ? - let (nrows, ncols) = self.col_piv_qr.data.shape(); + let (nrows, ncols) = self.col_piv_qr.shape_generic(); let mut res = OMatrix::identity_generic(nrows, ncols); if self.solve_mut(&mut res) { diff --git a/src/linalg/convolution.rs b/src/linalg/convolution.rs index 36cea3a0..21a32dbc 100644 --- a/src/linalg/convolution.rs +++ b/src/linalg/convolution.rs @@ -38,7 +38,7 @@ impl> Vector { .data .shape() .0 - .add(kernel.data.shape().0) + .add(kernel.shape_generic().0) .sub(Const::<1>); let mut conv = OVector::zeros_generic(result_len, Const::<1>); @@ -92,7 +92,7 @@ impl> Vector { .shape() .0 .add(Const::<1>) - .sub(kernel.data.shape().0); + .sub(kernel.shape_generic().0); let mut conv = OVector::zeros_generic(result_len, Const::<1>); for i in 0..(vec - ker + 1) { @@ -126,7 +126,7 @@ impl> Vector { panic!("convolve_same expects `self.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); } - let mut conv = OVector::zeros_generic(self.data.shape().0, Const::<1>); + let mut conv = OVector::zeros_generic(self.shape_generic().0, Const::<1>); for i in 0..vec { for j in 0..ker { diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs index c2816ff0..e7751af2 100644 --- a/src/linalg/exp.rs +++ b/src/linalg/exp.rs @@ -4,7 +4,6 @@ use crate::{ base::{ allocator::Allocator, dimension::{Const, Dim, DimMin, DimMinimum}, - storage::Storage, DefaultAllocator, }, convert, try_convert, ComplexField, OMatrix, RealField, @@ -47,7 +46,7 @@ where DefaultAllocator: Allocator + Allocator<(usize, usize), DimMinimum>, { fn new(a: OMatrix, use_exact_norm: bool) -> Self { - let (nrows, ncols) = a.data.shape(); + let (nrows, ncols) = a.shape_generic(); ExpmPadeHelper { use_exact_norm, ident: OMatrix::::identity_generic(nrows, ncols), @@ -348,7 +347,7 @@ where D: Dim, DefaultAllocator: Allocator + Allocator, { - let nrows = a.data.shape().0; + let nrows = a.shape_generic().0; let mut v = crate::OVector::::repeat_generic(nrows, Const::<1>, convert(1.0)); let m = a.transpose(); diff --git a/src/linalg/full_piv_lu.rs b/src/linalg/full_piv_lu.rs index f08af55c..20033c3c 100644 --- a/src/linalg/full_piv_lu.rs +++ b/src/linalg/full_piv_lu.rs @@ -53,7 +53,7 @@ where /// /// This effectively computes `P, L, U, Q` such that `P * matrix * Q = LU`. pub fn new(mut matrix: OMatrix) -> Self { - let (nrows, ncols) = matrix.data.shape(); + let (nrows, ncols) = matrix.shape_generic(); let min_nrows_ncols = nrows.min(ncols); let mut p = PermutationSequence::identity_generic(min_nrows_ncols); @@ -101,7 +101,7 @@ where where DefaultAllocator: Allocator>, { - let (nrows, ncols) = self.lu.data.shape(); + let (nrows, ncols) = self.lu.shape_generic(); let mut m = self.lu.columns_generic(0, nrows.min(ncols)).into_owned(); m.fill_upper_triangle(T::zero(), 1); m.fill_diagonal(T::one()); @@ -115,7 +115,7 @@ where where DefaultAllocator: Allocator, C>, { - let (nrows, ncols) = self.lu.data.shape(); + let (nrows, ncols) = self.lu.shape_generic(); self.lu.rows_generic(0, nrows.min(ncols)).upper_triangle() } @@ -222,7 +222,7 @@ where "FullPivLU inverse: unable to compute the inverse of a non-square matrix." ); - let (nrows, ncols) = self.lu.data.shape(); + let (nrows, ncols) = self.lu.shape_generic(); let mut res = OMatrix::identity_generic(nrows, ncols); if self.solve_mut(&mut res) { diff --git a/src/linalg/hessenberg.rs b/src/linalg/hessenberg.rs index 6b8ecfee..1e266b16 100644 --- a/src/linalg/hessenberg.rs +++ b/src/linalg/hessenberg.rs @@ -4,10 +4,11 @@ use serde::{Deserialize, Serialize}; use crate::allocator::Allocator; use crate::base::{DefaultAllocator, OMatrix, OVector}; use crate::dimension::{Const, DimDiff, DimSub, U1}; -use crate::storage::Storage; use simba::scalar::ComplexField; use crate::linalg::householder; +use crate::Matrix; +use std::mem::MaybeUninit; /// Hessenberg decomposition of a general matrix. #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] @@ -48,9 +49,7 @@ where { /// Computes the Hessenberg decomposition using householder reflections. pub fn new(hess: OMatrix) -> Self { - let mut work = unsafe { - crate::unimplemented_or_uninitialized_generic!(hess.data.shape().0, Const::<1>) - }; + let mut work = Matrix::zeros_generic(hess.shape_generic().0, Const::<1>); Self::new_with_workspace(hess, &mut work) } @@ -64,7 +63,7 @@ where "Cannot compute the hessenberg decomposition of a non-square matrix." ); - let dim = hess.data.shape().0; + let dim = hess.shape_generic().0; assert!( dim.value() != 0, @@ -76,18 +75,26 @@ where "Hessenberg: invalid workspace size." ); - let mut subdiag = unsafe { - crate::unimplemented_or_uninitialized_generic!(dim.sub(Const::<1>), Const::<1>) - }; - if dim.value() == 0 { - return Hessenberg { hess, subdiag }; + return Hessenberg { + hess, + subdiag: Matrix::zeros_generic(dim.sub(Const::<1>), Const::<1>), + }; } + let mut subdiag = Matrix::uninit(dim.sub(Const::<1>), Const::<1>); + for ite in 0..dim.value() - 1 { - householder::clear_column_unchecked(&mut hess, &mut subdiag[ite], ite, 1, Some(work)); + subdiag[ite] = MaybeUninit::new(householder::clear_column_unchecked( + &mut hess, + ite, + 1, + Some(work), + )); } + // Safety: subdiag is now fully initialized. + let subdiag = unsafe { subdiag.assume_init() }; Hessenberg { hess, subdiag } } diff --git a/src/linalg/householder.rs b/src/linalg/householder.rs index 9314ee45..6d20205d 100644 --- a/src/linalg/householder.rs +++ b/src/linalg/householder.rs @@ -3,7 +3,7 @@ use crate::allocator::Allocator; use crate::base::{DefaultAllocator, OMatrix, OVector, Unit, Vector}; use crate::dimension::Dim; -use crate::storage::{Storage, StorageMut}; +use crate::storage::StorageMut; use num::Zero; use simba::scalar::ComplexField; @@ -43,21 +43,23 @@ pub fn reflection_axis_mut>( /// Uses an householder reflection to zero out the `icol`-th column, starting with the `shift + 1`-th /// subdiagonal element. +/// +/// Returns the signed norm of the column. #[doc(hidden)] +#[must_use] pub fn clear_column_unchecked( matrix: &mut OMatrix, - diag_elt: &mut T, icol: usize, shift: usize, bilateral: Option<&mut OVector>, -) where +) -> T +where DefaultAllocator: Allocator + Allocator, { let (mut left, mut right) = matrix.columns_range_pair_mut(icol, icol + 1..); let mut axis = left.rows_range_mut(icol + shift..); let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis); - *diag_elt = reflection_norm; if not_zero { let refl = Reflection::new(Unit::new_unchecked(axis), T::zero()); @@ -67,19 +69,24 @@ pub fn clear_column_unchecked( } refl.reflect_with_sign(&mut right.rows_range_mut(icol + shift..), sign.conjugate()); } + + reflection_norm } /// Uses an householder reflection to zero out the `irow`-th row, ending before the `shift + 1`-th /// superdiagonal element. +/// +/// Returns the signed norm of the column. #[doc(hidden)] +#[must_use] pub fn clear_row_unchecked( matrix: &mut OMatrix, - diag_elt: &mut T, axis_packed: &mut OVector, work: &mut OVector, irow: usize, shift: usize, -) where +) -> T +where DefaultAllocator: Allocator + Allocator + Allocator, { let (mut top, mut bottom) = matrix.rows_range_pair_mut(irow, irow + 1..); @@ -88,7 +95,6 @@ pub fn clear_row_unchecked( let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis); axis.conjugate_mut(); // So that reflect_rows actually cancels the first row. - *diag_elt = reflection_norm; if not_zero { let refl = Reflection::new(Unit::new_unchecked(axis), T::zero()); @@ -102,6 +108,8 @@ pub fn clear_row_unchecked( } else { top.columns_range_mut(irow + shift..).tr_copy_from(&axis); } + + reflection_norm } /// Computes the orthogonal transformation described by the elementary reflector axii stored on @@ -113,7 +121,7 @@ where DefaultAllocator: Allocator, { assert!(m.is_square()); - let dim = m.data.shape().0; + let dim = m.shape_generic().0; // NOTE: we could build the identity matrix and call p_mult on it. // Instead we don't so that we take in account the matrix sparseness. diff --git a/src/linalg/lu.rs b/src/linalg/lu.rs index 36a00807..0e3be559 100644 --- a/src/linalg/lu.rs +++ b/src/linalg/lu.rs @@ -90,7 +90,7 @@ where { /// Computes the LU decomposition with partial (row) pivoting of `matrix`. pub fn new(mut matrix: OMatrix) -> Self { - let (nrows, ncols) = matrix.data.shape(); + let (nrows, ncols) = matrix.shape_generic(); let min_nrows_ncols = nrows.min(ncols); let mut p = PermutationSequence::identity_generic(min_nrows_ncols); @@ -132,7 +132,7 @@ where where DefaultAllocator: Allocator>, { - let (nrows, ncols) = self.lu.data.shape(); + let (nrows, ncols) = self.lu.shape_generic(); let mut m = self.lu.columns_generic(0, nrows.min(ncols)).into_owned(); m.fill_upper_triangle(T::zero(), 1); m.fill_diagonal(T::one()); @@ -149,7 +149,7 @@ where where DefaultAllocator: Reallocator>, { - let (nrows, ncols) = self.lu.data.shape(); + let (nrows, ncols) = self.lu.shape_generic(); let mut m = self.lu.resize_generic(nrows, nrows.min(ncols), T::zero()); m.fill_upper_triangle(T::zero(), 1); m.fill_diagonal(T::one()); @@ -162,7 +162,7 @@ where where DefaultAllocator: Reallocator>, { - let (nrows, ncols) = self.lu.data.shape(); + let (nrows, ncols) = self.lu.shape_generic(); let mut m = self.lu.resize_generic(nrows, nrows.min(ncols), T::zero()); m.fill_upper_triangle(T::zero(), 1); m.fill_diagonal(T::one()); @@ -176,7 +176,7 @@ where where DefaultAllocator: Allocator, C>, { - let (nrows, ncols) = self.lu.data.shape(); + let (nrows, ncols) = self.lu.shape_generic(); self.lu.rows_generic(0, nrows.min(ncols)).upper_triangle() } @@ -268,7 +268,7 @@ where "LU inverse: unable to compute the inverse of a non-square matrix." ); - let (nrows, ncols) = self.lu.data.shape(); + let (nrows, ncols) = self.lu.shape_generic(); let mut res = OMatrix::identity_generic(nrows, ncols); if self.try_inverse_to(&mut res) { Some(res) diff --git a/src/linalg/permutation_sequence.rs b/src/linalg/permutation_sequence.rs index ea868b5a..f4521988 100644 --- a/src/linalg/permutation_sequence.rs +++ b/src/linalg/permutation_sequence.rs @@ -69,11 +69,11 @@ where /// Creates a new sequence of D identity permutations. #[inline] pub fn identity_generic(dim: D) -> Self { - unsafe { - Self { - len: 0, - ipiv: crate::unimplemented_or_uninitialized_generic!(dim, Const::<1>), - } + Self { + len: 0, + // TODO: using a uninitialized matrix would save some computation, but + // that loos difficult to setup with MaybeUninit. + ipiv: Matrix::repeat_generic(dim, Const::<1>, (0, 0)), } } diff --git a/src/linalg/qr.rs b/src/linalg/qr.rs index 4bdbb364..e2f8e0c3 100644 --- a/src/linalg/qr.rs +++ b/src/linalg/qr.rs @@ -11,6 +11,7 @@ use simba::scalar::ComplexField; use crate::geometry::Reflection; use crate::linalg::householder; +use std::mem::MaybeUninit; /// The QR decomposition of a general matrix. #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] @@ -51,20 +52,25 @@ where { /// Computes the QR decomposition using householder reflections. pub fn new(mut matrix: OMatrix) -> Self { - let (nrows, ncols) = matrix.data.shape(); + let (nrows, ncols) = matrix.shape_generic(); let min_nrows_ncols = nrows.min(ncols); - let mut diag = - unsafe { crate::unimplemented_or_uninitialized_generic!(min_nrows_ncols, Const::<1>) }; - if min_nrows_ncols.value() == 0 { - return QR { qr: matrix, diag }; + return QR { + qr: matrix, + diag: Matrix::zeros_generic(min_nrows_ncols, Const::<1>), + }; } + let mut diag = Matrix::uninit(min_nrows_ncols, Const::<1>); + for i in 0..min_nrows_ncols.value() { - householder::clear_column_unchecked(&mut matrix, &mut diag[i], i, 0, None); + diag[i] = + MaybeUninit::new(householder::clear_column_unchecked(&mut matrix, i, 0, None)); } + // Safety: diag is now fully initialized. + let diag = unsafe { diag.assume_init() }; QR { qr: matrix, diag } } @@ -75,7 +81,7 @@ where where DefaultAllocator: Allocator, C>, { - let (nrows, ncols) = self.qr.data.shape(); + let (nrows, ncols) = self.qr.shape_generic(); let mut res = self.qr.rows_generic(0, nrows.min(ncols)).upper_triangle(); res.set_partial_diagonal(self.diag.iter().map(|e| T::from_real(e.modulus()))); res @@ -89,7 +95,7 @@ where where DefaultAllocator: Reallocator, C>, { - let (nrows, ncols) = self.qr.data.shape(); + let (nrows, ncols) = self.qr.shape_generic(); let mut res = self.qr.resize_generic(nrows.min(ncols), ncols, T::zero()); res.fill_lower_triangle(T::zero(), 1); res.set_partial_diagonal(self.diag.iter().map(|e| T::from_real(e.modulus()))); @@ -102,7 +108,7 @@ where where DefaultAllocator: Allocator>, { - let (nrows, ncols) = self.qr.data.shape(); + let (nrows, ncols) = self.qr.shape_generic(); // NOTE: we could build the identity matrix and call q_mul on it. // Instead we don't so that we take in account the matrix sparseness. @@ -254,7 +260,7 @@ where ); // TODO: is there a less naive method ? - let (nrows, ncols) = self.qr.data.shape(); + let (nrows, ncols) = self.qr.shape_generic(); let mut res = OMatrix::identity_generic(nrows, ncols); if self.solve_mut(&mut res) { diff --git a/src/linalg/schur.rs b/src/linalg/schur.rs index 3b650c52..953e9953 100644 --- a/src/linalg/schur.rs +++ b/src/linalg/schur.rs @@ -16,6 +16,8 @@ use crate::geometry::Reflection; use crate::linalg::givens::GivensRotation; use crate::linalg::householder; use crate::linalg::Hessenberg; +use crate::{Matrix, UninitVector}; +use std::mem::MaybeUninit; /// Schur decomposition of a square matrix. /// @@ -72,8 +74,7 @@ where /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm /// continues indefinitely until convergence. pub fn try_new(m: OMatrix, eps: T::RealField, max_niter: usize) -> Option { - let mut work = - unsafe { crate::unimplemented_or_uninitialized_generic!(m.data.shape().0, Const::<1>) }; + let mut work = Matrix::zeros_generic(m.shape_generic().0, Const::<1>); Self::do_decompose(m, &mut work, eps, max_niter, true) .map(|(q, t)| Schur { q: q.unwrap(), t }) @@ -91,7 +92,7 @@ where "Unable to compute the eigenvectors and eigenvalues of a non-square matrix." ); - let dim = m.data.shape().0; + let dim = m.shape_generic().0; // Specialization would make this easier. if dim.value() == 0 { @@ -294,7 +295,7 @@ where } /// Computes the complex eigenvalues of the decomposed matrix. - fn do_complex_eigenvalues(t: &OMatrix, out: &mut OVector, D>) + fn do_complex_eigenvalues(t: &OMatrix, out: &mut UninitVector, D>) where T: RealField, DefaultAllocator: Allocator, D>, @@ -306,7 +307,7 @@ where let n = m + 1; if t[(n, m)].is_zero() { - out[m] = NumComplex::new(t[(m, m)], T::zero()); + out[m] = MaybeUninit::new(NumComplex::new(t[(m, m)], T::zero())); m += 1; } else { // Solve the 2x2 eigenvalue subproblem. @@ -324,15 +325,15 @@ where let sqrt_discr = NumComplex::new(T::zero(), (-discr).sqrt()); let half_tra = (hnn + hmm) * crate::convert(0.5); - out[m] = NumComplex::new(half_tra, T::zero()) + sqrt_discr; - out[m + 1] = NumComplex::new(half_tra, T::zero()) - sqrt_discr; + out[m] = MaybeUninit::new(NumComplex::new(half_tra, T::zero()) + sqrt_discr); + out[m + 1] = MaybeUninit::new(NumComplex::new(half_tra, T::zero()) - sqrt_discr); m += 2; } } if m == dim - 1 { - out[m] = NumComplex::new(t[(m, m)], T::zero()); + out[m] = MaybeUninit::new(NumComplex::new(t[(m, m)], T::zero())); } } @@ -388,9 +389,7 @@ where /// Return `None` if some eigenvalues are complex. #[must_use] pub fn eigenvalues(&self) -> Option> { - let mut out = unsafe { - crate::unimplemented_or_uninitialized_generic!(self.t.data.shape().0, Const::<1>) - }; + let mut out = Matrix::zeros_generic(self.t.shape_generic().0, Const::<1>); if Self::do_eigenvalues(&self.t, &mut out) { Some(out) } else { @@ -405,11 +404,10 @@ where T: RealField, DefaultAllocator: Allocator, D>, { - let mut out = unsafe { - crate::unimplemented_or_uninitialized_generic!(self.t.data.shape().0, Const::<1>) - }; + let mut out = Matrix::uninit(self.t.shape_generic().0, Const::<1>); Self::do_complex_eigenvalues(&self.t, &mut out); - out + // Safety: out has been fully initialized by do_complex_eigenvalues. + unsafe { out.assume_init() } } } @@ -420,7 +418,7 @@ fn decompose_2x2( where DefaultAllocator: Allocator, { - let dim = m.data.shape().0; + let dim = m.shape_generic().0; let mut q = None; match compute_2x2_basis(&m.fixed_slice::<2, 2>(0, 0)) { Some(rot) => { @@ -519,9 +517,7 @@ where "Unable to compute eigenvalues of a non-square matrix." ); - let mut work = unsafe { - crate::unimplemented_or_uninitialized_generic!(self.data.shape().0, Const::<1>) - }; + let mut work = Matrix::zeros_generic(self.shape_generic().0, Const::<1>); // Special case for 2x2 matrices. if self.nrows() == 2 { @@ -547,6 +543,7 @@ where false, ) .unwrap(); + if Schur::do_eigenvalues(&schur.1, &mut work) { Some(work) } else { @@ -562,8 +559,8 @@ where T: RealField, DefaultAllocator: Allocator, D>, { - let dim = self.data.shape().0; - let mut work = unsafe { crate::unimplemented_or_uninitialized_generic!(dim, Const::<1>) }; + let dim = self.shape_generic().0; + let mut work = Matrix::zeros_generic(dim, Const::<1>); let schur = Schur::do_decompose( self.clone_owned(), @@ -573,8 +570,9 @@ where false, ) .unwrap(); - let mut eig = unsafe { crate::unimplemented_or_uninitialized_generic!(dim, Const::<1>) }; + let mut eig = Matrix::uninit(dim, Const::<1>); Schur::do_complex_eigenvalues(&schur.1, &mut eig); - eig + // Safety: eig has been fully initialized by do_complex_eigenvalues. + unsafe { eig.assume_init() } } } diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 241f00ce..0b50fc9b 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -111,7 +111,7 @@ where !matrix.is_empty(), "Cannot compute the SVD of an empty matrix." ); - let (nrows, ncols) = matrix.data.shape(); + let (nrows, ncols) = matrix.shape_generic(); let min_nrows_ncols = nrows.min(ncols); let dim = min_nrows_ncols.value(); diff --git a/src/linalg/symmetric_tridiagonal.rs b/src/linalg/symmetric_tridiagonal.rs index c7e87ba8..e071a916 100644 --- a/src/linalg/symmetric_tridiagonal.rs +++ b/src/linalg/symmetric_tridiagonal.rs @@ -4,10 +4,11 @@ use serde::{Deserialize, Serialize}; use crate::allocator::Allocator; use crate::base::{DefaultAllocator, OMatrix, OVector}; use crate::dimension::{Const, DimDiff, DimSub, U1}; -use crate::storage::Storage; use simba::scalar::ComplexField; use crate::linalg::householder; +use crate::Matrix; +use std::mem::MaybeUninit; /// Tridiagonalization of a symmetric matrix. #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] @@ -50,7 +51,7 @@ where /// /// Only the lower-triangular part (including the diagonal) of `m` is read. pub fn new(mut m: OMatrix) -> Self { - let dim = m.data.shape().0; + let dim = m.shape_generic().0; assert!( m.is_square(), @@ -61,19 +62,15 @@ where "Unable to compute the symmetric tridiagonal decomposition of an empty matrix." ); - let mut off_diagonal = unsafe { - crate::unimplemented_or_uninitialized_generic!(dim.sub(Const::<1>), Const::<1>) - }; - let mut p = unsafe { - crate::unimplemented_or_uninitialized_generic!(dim.sub(Const::<1>), Const::<1>) - }; + let mut off_diagonal = Matrix::uninit(dim.sub(Const::<1>), Const::<1>); + let mut p = Matrix::zeros_generic(dim.sub(Const::<1>), Const::<1>); for i in 0..dim.value() - 1 { let mut m = m.rows_range_mut(i + 1..); let (mut axis, mut m) = m.columns_range_pair_mut(i, i + 1..); let (norm, not_zero) = householder::reflection_axis_mut(&mut axis); - off_diagonal[i] = norm; + off_diagonal[i] = MaybeUninit::new(norm); if not_zero { let mut p = p.rows_range_mut(i..); @@ -87,6 +84,8 @@ where } } + // Safety: off_diagonal has been fully initialized. + let off_diagonal = unsafe { off_diagonal.assume_init() }; Self { tri: m, off_diagonal, diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs index 7b4a9cc9..546fa95a 100644 --- a/src/linalg/udu.rs +++ b/src/linalg/udu.rs @@ -4,7 +4,6 @@ use serde::{Deserialize, Serialize}; use crate::allocator::Allocator; use crate::base::{Const, DefaultAllocator, OMatrix, OVector}; use crate::dimension::Dim; -use crate::storage::Storage; use simba::scalar::RealField; /// UDU factorization. @@ -50,7 +49,7 @@ where /// Ref.: "Optimal control and estimation-Dover Publications", Robert F. Stengel, (1994) page 360 pub fn new(p: OMatrix) -> Option { let n = p.ncols(); - let n_dim = p.data.shape().1; + let n_dim = p.shape_generic().1; let mut d = OVector::zeros_generic(n_dim, Const::<1>); let mut u = OMatrix::zeros_generic(n_dim, n_dim); diff --git a/src/sparse/cs_matrix.rs b/src/sparse/cs_matrix.rs index 0fc3fed7..bb9f50a0 100644 --- a/src/sparse/cs_matrix.rs +++ b/src/sparse/cs_matrix.rs @@ -7,7 +7,7 @@ use std::slice; use crate::allocator::Allocator; use crate::sparse::cs_utils; -use crate::{Const, DefaultAllocator, Dim, Dynamic, OVector, Scalar, Vector, U1}; +use crate::{Const, DefaultAllocator, Dim, Dynamic, Matrix, OVector, Scalar, Vector, U1}; pub struct ColumnEntries<'a, T> { curr: usize, @@ -466,12 +466,12 @@ where { pub(crate) fn sort(&mut self) where + T: Zero, DefaultAllocator: Allocator, { // Size = R let nrows = self.data.shape().0; - let mut workspace = - unsafe { crate::unimplemented_or_uninitialized_generic!(nrows, Const::<1>) }; + let mut workspace = Matrix::zeros_generic(nrows, Const::<1>); self.sort_with_workspace(workspace.as_mut_slice()); } diff --git a/src/sparse/cs_matrix_cholesky.rs b/src/sparse/cs_matrix_cholesky.rs index 6d52d0a6..ff9ca023 100644 --- a/src/sparse/cs_matrix_cholesky.rs +++ b/src/sparse/cs_matrix_cholesky.rs @@ -3,7 +3,7 @@ use std::mem; use crate::allocator::Allocator; use crate::sparse::{CsMatrix, CsStorage, CsStorageIter, CsStorageIterMut, CsVecStorage}; -use crate::{Const, DefaultAllocator, Dim, OVector, RealField}; +use crate::{Const, DefaultAllocator, Dim, Matrix, OVector, RealField}; /// The cholesky decomposition of a column compressed sparse matrix. pub struct CsCholesky @@ -48,10 +48,8 @@ where let (l, u) = Self::nonzero_pattern(m); // Workspaces. - let work_x = - unsafe { crate::unimplemented_or_uninitialized_generic!(m.data.shape().0, Const::<1>) }; - let work_c = - unsafe { crate::unimplemented_or_uninitialized_generic!(m.data.shape().1, Const::<1>) }; + let work_x = Matrix::zeros_generic(m.data.shape().0, Const::<1>); + let work_c = Matrix::zeros_generic(m.data.shape().1, Const::<1>); let mut original_p = m.data.p.as_slice().to_vec(); original_p.push(m.data.i.len()); @@ -294,8 +292,7 @@ where let etree = Self::elimination_tree(m); let (nrows, ncols) = m.data.shape(); let mut rows = Vec::with_capacity(m.len()); - let mut cols = - unsafe { crate::unimplemented_or_uninitialized_generic!(m.data.shape().0, Const::<1>) }; + let mut cols = Matrix::zeros_generic(m.data.shape().0, Const::<1>); let mut marks = Vec::new(); // NOTE: the following will actually compute the non-zero pattern of diff --git a/src/sparse/cs_matrix_ops.rs b/src/sparse/cs_matrix_ops.rs index e03b12a5..419862a7 100644 --- a/src/sparse/cs_matrix_ops.rs +++ b/src/sparse/cs_matrix_ops.rs @@ -6,7 +6,7 @@ use crate::allocator::Allocator; use crate::constraint::{AreMultipliable, DimEq, ShapeConstraint}; use crate::sparse::{CsMatrix, CsStorage, CsStorageMut, CsVector}; use crate::storage::StorageMut; -use crate::{Const, DefaultAllocator, Dim, OVector, Scalar, Vector}; +use crate::{Const, DefaultAllocator, Dim, Matrix, OVector, Scalar, Vector}; impl> CsMatrix { fn scatter( @@ -219,7 +219,7 @@ where impl<'a, 'b, T, R1, R2, C1, C2, S1, S2> Add<&'b CsMatrix> for &'a CsMatrix where - T: Scalar + ClosedAdd + ClosedMul + One, + T: Scalar + ClosedAdd + ClosedMul + Zero + One, R1: Dim, C1: Dim, R2: Dim, @@ -242,8 +242,7 @@ where let mut res = CsMatrix::new_uninitialized_generic(nrows1, ncols2, self.len() + rhs.len()); let mut timestamps = OVector::zeros_generic(nrows1, Const::<1>); - let mut workspace = - unsafe { crate::unimplemented_or_uninitialized_generic!(nrows1, Const::<1>) }; + let mut workspace = Matrix::zeros_generic(nrows1, Const::<1>); let mut nz = 0; for j in 0..ncols2.value() { diff --git a/src/sparse/cs_matrix_solve.rs b/src/sparse/cs_matrix_solve.rs index 235fcef3..6136a0f8 100644 --- a/src/sparse/cs_matrix_solve.rs +++ b/src/sparse/cs_matrix_solve.rs @@ -152,8 +152,7 @@ impl> CsMatrix { self.lower_triangular_reach(b, &mut reach); // We sort the reach so the result matrix has sorted indices. reach.sort_unstable(); - let mut workspace = - unsafe { crate::unimplemented_or_uninitialized_generic!(b.data.shape().0, Const::<1>) }; + let mut workspace = Matrix::zeros_generic(b.data.shape().0, Const::<1>); for i in reach.iter().cloned() { workspace[i] = T::zero(); diff --git a/src/third_party/alga/alga_matrix.rs b/src/third_party/alga/alga_matrix.rs index e55ba49e..6a4cb982 100644 --- a/src/third_party/alga/alga_matrix.rs +++ b/src/third_party/alga/alga_matrix.rs @@ -15,8 +15,9 @@ use alga::linear::{ use crate::base::allocator::Allocator; use crate::base::dimension::{Dim, DimName}; -use crate::base::storage::{Storage, StorageMut}; -use crate::base::{DefaultAllocator, OMatrix, Scalar}; +use crate::base::storage::{RawStorage, RawStorageMut}; +use crate::base::{DefaultAllocator, Matrix, OMatrix, Scalar}; +use std::mem::MaybeUninit; /* * @@ -427,14 +428,14 @@ where { #[inline] fn meet_join(&self, other: &Self) -> (Self, Self) { - let shape = self.data.shape(); + let shape = self.shape_generic(); assert!( - shape == other.data.shape(), + shape == other.shape_generic(), "Matrix meet/join error: mismatched dimensions." ); - let mut mres = unsafe { crate::unimplemented_or_uninitialized_generic!(shape.0, shape.1) }; - let mut jres = unsafe { crate::unimplemented_or_uninitialized_generic!(shape.0, shape.1) }; + let mut mres = Matrix::uninit(shape.0, shape.1); + let mut jres = Matrix::uninit(shape.0, shape.1); for i in 0..shape.0.value() * shape.1.value() { unsafe { @@ -442,11 +443,12 @@ where .data .get_unchecked_linear(i) .meet_join(other.data.get_unchecked_linear(i)); - *mres.data.get_unchecked_linear_mut(i) = mj.0; - *jres.data.get_unchecked_linear_mut(i) = mj.1; + *mres.data.get_unchecked_linear_mut(i) = MaybeUninit::new(mj.0); + *jres.data.get_unchecked_linear_mut(i) = MaybeUninit::new(mj.1); } } - (mres, jres) + // Safety: both mres and jres are now completely initialized. + unsafe { (mres.assume_init(), jres.assume_init()) } } } diff --git a/src/third_party/glam/common/glam_matrix.rs b/src/third_party/glam/common/glam_matrix.rs index 77b68b5e..80f88054 100644 --- a/src/third_party/glam/common/glam_matrix.rs +++ b/src/third_party/glam/common/glam_matrix.rs @@ -2,7 +2,7 @@ use super::glam::{ BVec2, BVec3, BVec4, DMat2, DMat3, DMat4, DVec2, DVec3, DVec4, IVec2, IVec3, IVec4, Mat2, Mat3, Mat4, UVec2, UVec3, UVec4, Vec2, Vec3, Vec3A, Vec4, }; -use crate::storage::Storage; +use crate::storage::RawStorage; use crate::{Matrix, Matrix2, Matrix3, Matrix4, Vector, Vector2, Vector3, Vector4, U2, U3, U4}; macro_rules! impl_vec_conversion( @@ -16,7 +16,7 @@ macro_rules! impl_vec_conversion( impl From> for $Vec2 where - S: Storage<$N, U2>, + S: RawStorage<$N, U2>, { #[inline] fn from(e: Vector<$N, U2, S>) -> $Vec2 { @@ -33,7 +33,7 @@ macro_rules! impl_vec_conversion( impl From> for $Vec3 where - S: Storage<$N, U3>, + S: RawStorage<$N, U3>, { #[inline] fn from(e: Vector<$N, U3, S>) -> $Vec3 { @@ -50,7 +50,7 @@ macro_rules! impl_vec_conversion( impl From> for $Vec4 where - S: Storage<$N, U4>, + S: RawStorage<$N, U4>, { #[inline] fn from(e: Vector<$N, U4, S>) -> $Vec4 { @@ -75,7 +75,7 @@ impl From for Vector3 { impl From> for Vec3A where - S: Storage, + S: RawStorage, { #[inline] fn from(e: Vector) -> Vec3A { @@ -92,7 +92,7 @@ impl From for Matrix2 { impl From> for Mat2 where - S: Storage, + S: RawStorage, { #[inline] fn from(e: Matrix) -> Mat2 { @@ -112,7 +112,7 @@ impl From for Matrix3 { impl From> for Mat3 where - S: Storage, + S: RawStorage, { #[inline] fn from(e: Matrix) -> Mat3 { @@ -133,7 +133,7 @@ impl From for Matrix4 { impl From> for Mat4 where - S: Storage, + S: RawStorage, { #[inline] fn from(e: Matrix) -> Mat4 { @@ -155,7 +155,7 @@ impl From for Matrix2 { impl From> for DMat2 where - S: Storage, + S: RawStorage, { #[inline] fn from(e: Matrix) -> DMat2 { @@ -175,7 +175,7 @@ impl From for Matrix3 { impl From> for DMat3 where - S: Storage, + S: RawStorage, { #[inline] fn from(e: Matrix) -> DMat3 { @@ -196,7 +196,7 @@ impl From for Matrix4 { impl From> for DMat4 where - S: Storage, + S: RawStorage, { #[inline] fn from(e: Matrix) -> DMat4 { diff --git a/src/third_party/mint/mint_matrix.rs b/src/third_party/mint/mint_matrix.rs index 1e0a4d54..ce45fcda 100644 --- a/src/third_party/mint/mint_matrix.rs +++ b/src/third_party/mint/mint_matrix.rs @@ -1,10 +1,10 @@ use std::convert::{AsMut, AsRef, From, Into}; -use std::mem; +use std::mem::{self, MaybeUninit}; use std::ptr; use crate::base::allocator::Allocator; -use crate::base::dimension::{U1, U2, U3, U4}; -use crate::base::storage::{ContiguousStorage, ContiguousStorageMut, Storage, StorageMut}; +use crate::base::dimension::{Const, DimName, U1, U2, U3, U4}; +use crate::base::storage::{IsContiguous, RawStorage, RawStorageMut}; use crate::base::{DefaultAllocator, Matrix, OMatrix, Scalar}; macro_rules! impl_from_into_mint_1D( @@ -15,9 +15,12 @@ macro_rules! impl_from_into_mint_1D( #[inline] fn from(v: mint::$VT) -> Self { unsafe { - let mut res = Self::new_uninitialized(); - ptr::copy_nonoverlapping(&v.x, (*res.as_mut_ptr()).data.ptr_mut(), $SZ); - + let mut res = Matrix::uninit(<$NRows>::name(), Const::<1>); + // Copy the data. + ptr::copy_nonoverlapping(&v.x, res.data.ptr_mut() as *mut T, $SZ); + // Prevent from being dropped the originals we just copied. + mem::forget(v); + // The result is now fully initialized. res.assume_init() } } @@ -25,22 +28,28 @@ macro_rules! impl_from_into_mint_1D( impl Into> for Matrix where T: Scalar, - S: ContiguousStorage { + S: RawStorage + IsContiguous { #[inline] fn into(self) -> mint::$VT { + // SAFETY: this is OK thanks to the IsContiguous bound. unsafe { - let mut res: mint::$VT = mem::MaybeUninit::uninit().assume_init(); - ptr::copy_nonoverlapping(self.data.ptr(), &mut res.x, $SZ); - res + let mut res: MaybeUninit> = MaybeUninit::uninit(); + // Copy the data. + ptr::copy_nonoverlapping(self.data.ptr(), res.as_mut_ptr() as *mut T, $SZ); + // Prevent from being dropped the originals we just copied. + mem::forget(self); + // The result is now fully initialized. + res.assume_init() } } } impl AsRef> for Matrix where T: Scalar, - S: ContiguousStorage { + S: RawStorage + IsContiguous { #[inline] fn as_ref(&self) -> &mint::$VT { + // SAFETY: this is OK thanks to the IsContiguous bound. unsafe { mem::transmute(self.data.ptr()) } @@ -49,9 +58,10 @@ macro_rules! impl_from_into_mint_1D( impl AsMut> for Matrix where T: Scalar, - S: ContiguousStorageMut { + S: RawStorageMut + IsContiguous { #[inline] fn as_mut(&mut self) -> &mut mint::$VT { + // SAFETY: this is OK thanks to the IsContiguous bound. unsafe { mem::transmute(self.data.ptr_mut()) } @@ -75,13 +85,15 @@ macro_rules! impl_from_into_mint_2D( #[inline] fn from(m: mint::$MV) -> Self { unsafe { - let mut res = Self::new_uninitialized(); - let mut ptr = (*res.as_mut_ptr()).data.ptr_mut(); + let mut res = Matrix::uninit(<$NRows>::name(), <$NCols>::name()); + let mut ptr = res.data.ptr_mut(); $( - ptr::copy_nonoverlapping(&m.$component.x, ptr, $SZRows); + ptr::copy_nonoverlapping(&m.$component.x, ptr as *mut T, $SZRows); ptr = ptr.offset($SZRows); )* - let _ = ptr; + let _ = ptr; // Just to avoid some unused assignment warnings. + // Forget the original data to avoid double-free. + mem::forget(m); res.assume_init() } } @@ -93,14 +105,16 @@ macro_rules! impl_from_into_mint_2D( #[inline] fn into(self) -> mint::$MV { unsafe { - let mut res: mint::$MV = mem::MaybeUninit::uninit().assume_init(); + let mut res: MaybeUninit> = MaybeUninit::uninit(); let mut ptr = self.data.ptr(); $( - ptr::copy_nonoverlapping(ptr, &mut res.$component.x, $SZRows); + ptr::copy_nonoverlapping(ptr, ptr::addr_of_mut!((*res.as_mut_ptr()).$component) as *mut T, $SZRows); ptr = ptr.offset($SZRows); )* let _ = ptr; - res + // Forget the original data to avoid double-free. + mem::forget(self); + res.assume_init() } } } diff --git a/src/third_party/mint/mint_point.rs b/src/third_party/mint/mint_point.rs index fbce1c88..45f85e3c 100644 --- a/src/third_party/mint/mint_point.rs +++ b/src/third_party/mint/mint_point.rs @@ -1,4 +1,4 @@ -use crate::base::storage::{Storage, StorageMut}; +use crate::base::storage::{RawStorage, RawStorageMut}; use crate::{OVector, Point, Scalar}; use std::convert::{AsMut, AsRef}; diff --git a/tests/core/matrix.rs b/tests/core/matrix.rs index eaa252db..4a35fb20 100644 --- a/tests/core/matrix.rs +++ b/tests/core/matrix.rs @@ -447,7 +447,7 @@ fn apply() { 1.0, 2.0, 3.0, 4.0, 6.0, 7.0, 8.0, 9.0, 10.0, 9.0, 8.0, 7.0, 6.0, 4.0, 3.0, 2.0, ); - a.apply(|e| e.round()); + a.apply(|e| *e = e.round()); assert_eq!(a, expected); }