From df9b6f5f646e90eb6300e5b08a39e253e5474e88 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Violeta=20Hern=C3=A1ndez?= Date: Thu, 15 Jul 2021 23:56:58 -0500 Subject: [PATCH] blas.rs works now! --- src/base/allocator.rs | 12 +- src/base/blas.rs | 407 ++++++++++++++++------------- src/base/construction.rs | 97 +++++-- src/base/conversion.rs | 14 +- src/base/default_allocator.rs | 12 +- src/base/edition.rs | 70 +++-- src/base/matrix.rs | 232 ++++++++++++---- src/base/matrix_slice.rs | 64 ++++- src/base/ops.rs | 176 ++++++------- src/base/statistics.rs | 14 +- src/linalg/permutation_sequence.rs | 12 +- 11 files changed, 695 insertions(+), 415 deletions(-) diff --git a/src/base/allocator.rs b/src/base/allocator.rs index 77c9b528..92a38300 100644 --- a/src/base/allocator.rs +++ b/src/base/allocator.rs @@ -49,7 +49,7 @@ pub trait Allocator: /// A matrix reallocator. Changes the size of the memory buffer that initially contains (RFrom × /// CFrom) elements to a smaller or larger size (RTo, CTo). pub trait Reallocator: - InnerAllocator + InnerAllocator + Allocator + Allocator { /// Reallocates a buffer of shape `(RTo, CTo)`, possibly reusing a previously allocated buffer /// `buf`. Data stored by `buf` are linearly copied to the output: @@ -75,7 +75,7 @@ pub type SameShapeC = >:: // TODO: Bad name. /// Restricts the given number of rows and columns to be respectively the same. pub trait SameShapeAllocator: - InnerAllocator + InnerAllocator, SameShapeC> + Allocator + Allocator, SameShapeC> where ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { @@ -85,7 +85,7 @@ impl SameShapeAllocator + InnerAllocator, SameShapeC>, + Allocator + Allocator, SameShapeC>, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { } @@ -93,8 +93,8 @@ where // XXX: Bad name. /// Restricts the given number of rows to be equal. pub trait SameShapeVectorAllocator: - InnerAllocator - + InnerAllocator> + Allocator + + Allocator> + SameShapeAllocator where ShapeConstraint: SameNumberOfRows, @@ -103,7 +103,7 @@ where impl SameShapeVectorAllocator for DefaultAllocator where - DefaultAllocator: InnerAllocator + InnerAllocator>, + DefaultAllocator: Allocator + Allocator>, ShapeConstraint: SameNumberOfRows, { } diff --git a/src/base/blas.rs b/src/base/blas.rs index 3b8ac951..2ef0dff7 100644 --- a/src/base/blas.rs +++ b/src/base/blas.rs @@ -1,4 +1,12 @@ -use crate::{OVector, SimdComplexField}; +//! Implements a subset of the Basic Linear Algebra Subprograms (BLAS), a +//! standard and highly optimized set of basic vector and matrix operations. +//! +//! To avoid unsoundness due to mishandling of uninitialized data, we divide our +//! methods into two groups: those that take in a `&mut` to a matrix, and those +//! that return an owned matrix that would otherwise result from setting a +//! parameter to zero in the other methods. + +use crate::{OMatrix, OVector, SimdComplexField}; #[cfg(feature = "std")] use matrixmultiply; use num::{One, Zero}; @@ -279,72 +287,16 @@ 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(); - } - } -} - -fn array_axc_uninit( - y: &mut [MaybeUninit], - 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) = MaybeUninit::new( - a.inlined_clone() - * x.get_unchecked(i * stride2).inlined_clone() - * c.inlined_clone(), - ); - } - } -} - /// # BLAS functions impl Vector where T: Scalar + Zero + ClosedAdd + ClosedMul, S: StorageMut, { - /// Computes `self = a * x * c + b * self`. + /// Computes `self = a * x * c + b * self`, where `a`, `b`, `c` are scalars, + /// and `x` is a vector of the same size as `self`. + /// + /// For commutative scalars, this is equivalent to an [`axpy`] call. /// /// If `b` is zero, `self` is never read from. /// @@ -376,9 +328,24 @@ where let x = x.data.as_slice_unchecked(); if !b.is_zero() { - array_axcpy(y, a, x, c, b, rstride1, rstride2, x.len()); + for i in 0..x.len() { + unsafe { + let y = y.get_unchecked_mut(i * rstride1); + *y = a.inlined_clone() + * x.get_unchecked(i * rstride2).inlined_clone() + * c.inlined_clone() + + b.inlined_clone() * y.inlined_clone(); + } + } } else { - array_axc(y, a, x, c, rstride1, rstride2, x.len()); + for i in 0..x.len() { + unsafe { + let y = y.get_unchecked_mut(i * rstride1); + *y = a.inlined_clone() + * x.get_unchecked(i * rstride2).inlined_clone() + * c.inlined_clone(); + } + } } } } @@ -746,49 +713,55 @@ where } } -impl OVector, D> +impl Vector, D, S> where T: Scalar + Zero + ClosedAdd + ClosedMul, - DefaultAllocator: Allocator, + S: StorageMut, D>, { - pub fn axc(&mut self, a: T, x: &Vector, c: T) -> OVector + /// Computes `alpha * a * x`, where `a` is a matrix, `x` a vector, and + /// `alpha` is a scalar. + /// + /// # Safety + /// `self` must be completely uninitialized, or data leaks will occur. After + /// this method is called, all entries in `self` will be initialized. + pub fn axc(&mut self, a: T, x: &Vector, c: T) where - SB: Storage, + S2: 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(); - array_axc_uninit(y, a, x, c, rstride1, rstride2, x.len()); - self.assume_init() + for i in 0..y.len() { + *y.get_unchecked_mut(i * rstride1) = MaybeUninit::new( + a.inlined_clone() + * x.get_unchecked(i * rstride2).inlined_clone() + * c.inlined_clone(), + ); + } } } - /// Computes `self = alpha * a * x, where `a` is a matrix, `x` a vector, and + /// Computes `alpha * a * x`, where `a` is a matrix, `x` a vector, and /// `alpha` is a scalar. /// - /// By the time this method returns, `self` will have been initialized. + /// Initializes `self`. #[inline] - pub fn gemv_uninit( - mut self, + pub fn gemv_z( + &mut self, alpha: T, a: &Matrix, x: &Vector, - beta: T, - ) -> OVector - where + ) where T: One, SB: Storage, SC: Storage, ShapeConstraint: DimEq + AreMultipliable, + // DefaultAllocator: Allocator, { let dim1 = self.nrows(); let (nrows2, ncols2) = a.shape(); @@ -801,22 +774,169 @@ where if ncols2 == 0 { self.fill_fn(|| MaybeUninit::new(T::zero())); - return self.assume_init(); + return; } // TODO: avoid bound checks. let col2 = a.column(0); let val = unsafe { x.vget_unchecked(0).inlined_clone() }; - let res = self.axc(alpha.inlined_clone(), &col2, val); + self.axc(alpha.inlined_clone(), &col2, val); - for j in 1..ncols2 { - let col2 = a.column(j); - let val = unsafe { x.vget_unchecked(j).inlined_clone() }; + // Safety: axc initializes self. + unsafe { + let mut init = self.assume_init_mut(); - res.axcpy(alpha.inlined_clone(), &col2, val, T::one()); + for j in 1..ncols2 { + let col2 = a.column(j); + let val = unsafe { x.vget_unchecked(j).inlined_clone() }; + init.axcpy(alpha.inlined_clone(), &col2, val, T::one()); + } + } + } +} + +impl OMatrix +where + T: Scalar + Zero + One + ClosedAdd + ClosedMul, + DefaultAllocator: Allocator, +{ + /// Computes `alpha * a * b`, where `a` and `b` are matrices, and `alpha` is + /// a scalar. + /// + /// # Examples: + /// + /// ``` + /// # #[macro_use] extern crate approx; + /// # use nalgebra::{Matrix2x3, Matrix3x4, Matrix2x4}; + /// let mut mat1 = Matrix2x4::identity(); + /// let mat2 = Matrix2x3::new(1.0, 2.0, 3.0, + /// 4.0, 5.0, 6.0); + /// let mat3 = Matrix3x4::new(0.1, 0.2, 0.3, 0.4, + /// 0.5, 0.6, 0.7, 0.8, + /// 0.9, 1.0, 1.1, 1.2); + /// let expected = mat2 * mat3 * 10.0 + mat1 * 5.0; + /// + /// mat1.gemm(10.0, &mat2, &mat3, 5.0); + /// assert_relative_eq!(mat1, expected); + /// ``` + #[inline] + pub fn gemm_z( + alpha: T, + a: &Matrix, + b: &Matrix, + ) -> Self + where + SB: Storage, + SC: Storage, + ShapeConstraint: SameNumberOfRows + + SameNumberOfColumns + + AreMultipliable, + { + let (nrows1, ncols1) = a.shape(); + let (nrows2, ncols2) = b.shape(); + + assert_eq!( + ncols1, nrows2, + "gemm: dimensions mismatch for multiplication." + ); + + let mut res = + Matrix::new_uninitialized_generic(R1::from_usize(nrows1), C1::from_usize(ncols2)); + + #[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. + + // Threshold determined empirically. + const SMALL_DIM: usize = 5; + + if nrows1 > SMALL_DIM + && ncols1 > SMALL_DIM + && nrows2 > SMALL_DIM + && ncols2 > SMALL_DIM + { + // 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 ncols1 == 0 { + // NOTE: we can't just always multiply by beta + // because we documented the guaranty that `self` is + // never read if `beta` is zero. + + // Safety: this buffer is empty. + return res.assume_init(); + } + + let (rsa, csa) = a.strides(); + let (rsb, csb) = b.strides(); + let (rsc, csc) = res.strides(); + + if T::is::() { + unsafe { + matrixmultiply::sgemm( + nrows1, + ncols1, + ncols2, + 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, + 0.0, + res.data.ptr_mut() as *mut f32, + rsc as isize, + csc as isize, + ); + + return res.assume_init(); + } + } else if T::is::() { + unsafe { + matrixmultiply::dgemm( + nrows1, + ncols1, + ncols2, + 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, + 0.0, + res.data.ptr_mut() as *mut f64, + rsc as isize, + csc as isize, + ); + + return res.assume_init(); + } + } + } + } } - res + for j1 in 0..ncols1 { + // TODO: avoid bound checks. + res.column_mut(j1) + .gemv_z(alpha.inlined_clone(), a, &b.column(j1)); + } + + unsafe { res.assume_init() } } } @@ -1372,49 +1492,6 @@ where /// /// mat.quadform_tr_with_workspace(&mut workspace, 10.0, &lhs, &mid, 5.0); /// assert_relative_eq!(mat, expected); - pub fn quadform_tr_with_workspace( - &mut self, - work: &mut OVector, D2>, - alpha: T, - lhs: &Matrix, - mid: &SquareMatrix, - beta: T, - ) where - S3: Storage, - S4: Storage, - ShapeConstraint: DimEq + DimEq + DimEq + DimEq, - DefaultAllocator: Allocator, - { - let work = work.gemv_uninit(T::one(), lhs, &mid.column(0), T::zero()); - self.ger(alpha.inlined_clone(), &work, &lhs.column(0), beta); - - for j in 1..mid.ncols() { - work.gemv(T::one(), lhs, &mid.column(j), T::zero()); - self.ger(alpha.inlined_clone(), &work, &lhs.column(j), T::one()); - } - } - - /// Computes the quadratic form `self = alpha * lhs * mid * lhs.transpose() + beta * self`. - /// - /// This allocates a workspace vector of dimension D1 for intermediate results. - /// If `D1` is a type-level integer, then the allocation is performed on the stack. - /// Use `.quadform_tr_with_workspace(...)` instead to avoid allocations. - /// - /// # Examples: - /// - /// ``` - /// # #[macro_use] extern crate approx; - /// # use nalgebra::{Matrix2, Matrix3, Matrix2x3, Vector2}; - /// let mut mat = Matrix2::identity(); - /// let lhs = Matrix2x3::new(1.0, 2.0, 3.0, - /// 4.0, 5.0, 6.0); - /// let mid = Matrix3::new(0.1, 0.2, 0.3, - /// 0.5, 0.6, 0.7, - /// 0.9, 1.0, 1.1); - /// let expected = lhs * mid * lhs.transpose() * 10.0 + mat * 5.0; - /// - /// mat.quadform_tr(10.0, &lhs, &mid, 5.0); - /// assert_relative_eq!(mat, expected); pub fn quadform_tr( &mut self, alpha: T, @@ -1424,11 +1501,19 @@ where ) where S3: Storage, S4: Storage, - ShapeConstraint: DimEq + DimEq + DimEq, - DefaultAllocator: Allocator, + ShapeConstraint: DimEq + DimEq, + DefaultAllocator: Allocator, { - let mut work = Matrix::new_uninitialized_generic(self.data.shape().0, Const::<1>); - self.quadform_tr_with_workspace(&mut work, alpha, lhs, mid, beta) + let work = Matrix::new_uninitialized_generic(R3::from_usize(self.shape().0), Const::<1>); + work.gemv_z(T::one(), lhs, &mid.column(0)); + let work = unsafe { work.assume_init() }; + + self.ger(alpha.inlined_clone(), &work, &lhs.column(0), beta); + + for j in 1..mid.ncols() { + work.gemv(T::one(), lhs, &mid.column(j), T::zero()); + self.ger(alpha.inlined_clone(), &work, &lhs.column(j), T::one()); + } } /// Computes the quadratic form `self = alpha * rhs.transpose() * mid * rhs + beta * self`. @@ -1454,11 +1539,10 @@ where /// let mut workspace = DVector::new_random(3); /// let expected = rhs.transpose() * &mid * &rhs * 10.0 + &mat * 5.0; /// - /// mat.quadform_with_workspace(&mut workspace, 10.0, &mid, &rhs, 5.0); + /// mat.quadform(&mut workspace, 10.0, &mid, &rhs, 5.0); /// assert_relative_eq!(mat, expected); - pub fn quadform_with_workspace( + pub fn quadform( &mut self, - work: &mut OVector, D2>, alpha: T, mid: &SquareMatrix, rhs: &Matrix, @@ -1466,54 +1550,21 @@ where ) where S3: Storage, S4: Storage, - ShapeConstraint: - DimEq + DimEq + DimEq + AreMultipliable, - DefaultAllocator: Allocator, + ShapeConstraint: DimEq + DimEq + DimEq, + DefaultAllocator: Allocator, { - let work = work.gemv_uninit(T::one(), mid, &rhs.column(0), T::zero()); + // TODO: figure out why type inference wasn't doing its job. + let work = Matrix::new_uninitialized_generic(D3::from_usize(self.shape().0), Const::<1>); + work.gemv_z::(T::one(), mid, &rhs.column(0)); + let work = unsafe { work.assume_init() }; + self.column_mut(0) .gemv_tr(alpha.inlined_clone(), rhs, &work, beta.inlined_clone()); for j in 1..rhs.ncols() { - work.gemv(T::one(), mid, &rhs.column(j), T::zero()); + work.gemv::(T::one(), mid, &rhs.column(j), T::zero()); self.column_mut(j) .gemv_tr(alpha.inlined_clone(), rhs, &work, beta.inlined_clone()); } } - - /// Computes the quadratic form `self = alpha * rhs.transpose() * mid * rhs + beta * self`. - /// - /// This allocates a workspace vector of dimension D2 for intermediate results. - /// If `D2` is a type-level integer, then the allocation is performed on the stack. - /// Use `.quadform_with_workspace(...)` instead to avoid allocations. - /// - /// ``` - /// # #[macro_use] extern crate approx; - /// # use nalgebra::{Matrix2, Matrix3x2, Matrix3}; - /// let mut mat = Matrix2::identity(); - /// let rhs = Matrix3x2::new(1.0, 2.0, - /// 3.0, 4.0, - /// 5.0, 6.0); - /// let mid = Matrix3::new(0.1, 0.2, 0.3, - /// 0.5, 0.6, 0.7, - /// 0.9, 1.0, 1.1); - /// let expected = rhs.transpose() * mid * rhs * 10.0 + mat * 5.0; - /// - /// mat.quadform(10.0, &mid, &rhs, 5.0); - /// assert_relative_eq!(mat, expected); - pub fn quadform( - &mut self, - alpha: T, - mid: &SquareMatrix, - rhs: &Matrix, - beta: T, - ) where - S2: Storage, - S3: Storage, - ShapeConstraint: DimEq + DimEq + AreMultipliable, - DefaultAllocator: Allocator, - { - let mut work = Matrix::new_uninitialized_generic(mid.data.shape().0, Const::<1>); - self.quadform_with_workspace(&mut work, alpha, mid, rhs, beta) - } } diff --git a/src/base/construction.rs b/src/base/construction.rs index c040a9dc..f0709917 100644 --- a/src/base/construction.rs +++ b/src/base/construction.rs @@ -18,15 +18,14 @@ use typenum::{self, Cmp, Greater}; use simba::scalar::{ClosedAdd, ClosedMul}; -use crate::{base::allocator::Allocator}; +use crate::base::allocator::{Allocator, InnerAllocator}; use crate::base::dimension::{Dim, DimName, Dynamic, ToTypenum}; use crate::base::storage::Storage; use crate::base::{ ArrayStorage, Const, DefaultAllocator, Matrix, OMatrix, OVector, Scalar, Unit, Vector, }; -/// 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 +/// OBJECTIVE: GET RID OF THIS! #[macro_export] macro_rules! unimplemented_or_uninitialized_generic { ($nrows:expr, $ncols:expr) => {{ @@ -99,7 +98,7 @@ where "Matrix init. error: the slice did not contain the right number of elements." ); - let mut res = Matrix::new_uninitialized_generic(nrows, ncols); + let mut res = Self::new_uninitialized_generic(nrows, ncols); let mut iter = slice.iter(); for i in 0..nrows.value() { @@ -117,7 +116,10 @@ where /// Creates a matrix with its elements filled with the components provided by a slice. The /// components must have the same layout as the matrix data storage (i.e. column-major). #[inline] - pub fn from_column_slice_generic(nrows: R, ncols: C, slice: &[T]) -> Self where T:Clone{ + pub fn from_column_slice_generic(nrows: R, ncols: C, slice: &[T]) -> Self + where + T: Clone, + { Self::from_iterator_generic(nrows, ncols, slice.iter().cloned()) } @@ -128,7 +130,7 @@ where where F: FnMut(usize, usize) -> T, { - let mut res = Matrix::new_uninitialized_generic(nrows, ncols); + let mut res = Self::new_uninitialized_generic(nrows, ncols); for j in 0..ncols.value() { for i in 0..nrows.value() { @@ -139,7 +141,7 @@ where } // Safety: all entries have been initialized. - unsafe { res.assume_init()} + unsafe { res.assume_init() } } /// Creates a new identity matrix. @@ -352,7 +354,7 @@ where #[inline] pub fn from_diagonal>(diag: &Vector) -> Self where - T: Zero+Scalar, + T: Zero + Scalar, { let (dim, _) = diag.data.shape(); let mut res = Self::zeros_generic(dim, dim); @@ -374,12 +376,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),*) -> MaybeUninit { - Self::new_uninitialized_generic($($gargs),*) - } - /// Creates a matrix or vector with all its elements set to `elem`. /// /// # Example @@ -518,8 +514,7 @@ macro_rules! impl_constructors( /// dm[(1, 0)] == 3 && dm[(1, 1)] == 4 && dm[(1, 2)] == 5); /// ``` #[inline] - pub fn from_fn($($args: usize,)* f: F) -> Self - where F: FnMut(usize, usize) -> T { + pub fn from_fn T>($($args: usize,)* f: F) -> Self { Self::from_fn_generic($($gargs, )* f) } @@ -543,7 +538,9 @@ macro_rules! impl_constructors( /// ``` #[inline] pub fn identity($($args: usize,)*) -> Self - where T: Zero + One { + where + T: Zero + One + Scalar + { Self::identity_generic($($gargs),* ) } @@ -566,7 +563,9 @@ macro_rules! impl_constructors( /// ``` #[inline] pub fn from_diagonal_element($($args: usize,)* elt: T) -> Self - where T: Zero + One { + where + T: Zero + One + Scalar + { Self::from_diagonal_element_generic($($gargs, )* elt) } @@ -593,7 +592,9 @@ macro_rules! impl_constructors( /// ``` #[inline] pub fn from_partial_diagonal($($args: usize,)* elts: &[T]) -> Self - where T: Zero { + where + T: Zero + Scalar + { Self::from_partial_diagonal_generic($($gargs, )* elts) } @@ -612,7 +613,9 @@ macro_rules! impl_constructors( #[inline] #[cfg(feature = "rand")] pub fn new_random($($args: usize),*) -> Self - where Standard: Distribution { + where + Standard: Distribution + { Self::new_random_generic($($gargs),*) } } @@ -630,6 +633,17 @@ where ); // Arguments for non-generic constructors. } +impl OMatrix, R, C> +where + DefaultAllocator: Allocator, +{ + /// Creates a new uninitialized matrix or vector. + #[inline] + pub fn new_uninitialized() -> Self { + Self::new_uninitialized_generic(R::name(), C::name()) + } +} + /// # Constructors of matrices with a dynamic number of columns impl OMatrix where @@ -641,6 +655,17 @@ where ncols); } +impl OMatrix, R, Dynamic> +where + DefaultAllocator: Allocator, +{ + /// Creates a new uninitialized matrix or vector. + #[inline] + pub fn new_uninitialized(ncols: usize) -> Self { + Self::new_uninitialized_generic(R::name(), Dynamic::new(ncols)) + } +} + /// # Constructors of dynamic vectors and matrices with a dynamic number of rows impl OMatrix where @@ -652,6 +677,17 @@ where nrows); } +impl OMatrix, Dynamic, C> +where + DefaultAllocator: Allocator, +{ + /// Creates a new uninitialized matrix or vector. + #[inline] + pub fn new_uninitialized(nrows: usize) -> Self { + Self::new_uninitialized_generic(Dynamic::new(nrows), C::name()) + } +} + /// # Constructors of fully dynamic matrices impl OMatrix where @@ -663,6 +699,17 @@ where nrows, ncols); } +impl OMatrix, Dynamic, Dynamic> +where + DefaultAllocator: Allocator, +{ + /// Creates a new uninitialized matrix or vector. + #[inline] + pub fn new_uninitialized(nrows: usize, ncols: usize) -> Self { + Self::new_uninitialized_generic(Dynamic::new(nrows), Dynamic::new(ncols)) + } +} + /* * * Constructors that don't necessarily require all dimensions @@ -701,7 +748,10 @@ macro_rules! impl_constructors_from_data( /// dm[(1, 0)] == 3 && dm[(1, 1)] == 4 && dm[(1, 2)] == 5); /// ``` #[inline] - pub fn from_row_slice($($args: usize,)* $data: &[T]) -> Self { + pub fn from_row_slice($($args: usize,)* $data: &[T]) -> Self + where + T: Clone + { Self::from_row_slice_generic($($gargs, )* $data) } @@ -728,7 +778,10 @@ macro_rules! impl_constructors_from_data( /// dm[(1, 0)] == 1 && dm[(1, 1)] == 3 && dm[(1, 2)] == 5); /// ``` #[inline] - pub fn from_column_slice($($args: usize,)* $data: &[T]) -> Self { + pub fn from_column_slice($($args: usize,)* $data: &[T]) -> Self + where + T: Clone + { Self::from_column_slice_generic($($gargs, )* $data) } diff --git a/src/base/conversion.rs b/src/base/conversion.rs index 97194a13..1efb9a91 100644 --- a/src/base/conversion.rs +++ b/src/base/conversion.rs @@ -27,14 +27,10 @@ use crate::constraint::DimEq; use crate::{IsNotStaticOne, RowSVector, SMatrix, SVector}; // TODO: too bad this won't work for slice conversions. -impl SubsetOf> for OMatrix +impl SubsetOf> + for OMatrix where - R1: Dim, - C1: Dim, - R2: Dim, - C2: Dim, - T1: Scalar, - T2: Scalar + SupersetOf, + T2: SupersetOf, DefaultAllocator: Allocator + Allocator + SameShapeAllocator, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, @@ -45,7 +41,7 @@ where let nrows2 = R2::from_usize(nrows); let ncols2 = C2::from_usize(ncols); - let mut res = OMatrix::::new_uninitialized_generic(nrows2, ncols2); + let mut res = Matrix::new_uninitialized_generic(nrows2, ncols2); for i in 0..nrows { for j in 0..ncols { @@ -57,7 +53,7 @@ where } // Safety: all entries have been initialized. - unsafe { Matrix::assume_init(res) } + unsafe { res.assume_init() } } #[inline] diff --git a/src/base/default_allocator.rs b/src/base/default_allocator.rs index 7ee425ff..b9cb793c 100644 --- a/src/base/default_allocator.rs +++ b/src/base/default_allocator.rs @@ -77,9 +77,13 @@ impl Allocator, Const> for Def unsafe fn assume_init( uninit: , Const, Const>>::Buffer, ) -> Owned, Const> { - // Safety: MaybeUninit has the same alignment and layout as T, and by - // extension so do arrays based on these. - mem::transmute(uninit) + // SAFETY: + // * The caller guarantees that all elements of the array are initialized + // * `MaybeUninit` and T are guaranteed to have the same layout + // * MaybeUnint does not drop, so there are no double-frees + // * `ArrayStorage` is transparent. + // And thus the conversion is safe + ArrayStorage((&uninit as *const _ as *const [_; C]).read()) } } @@ -205,7 +209,7 @@ where ); // Safety: TODO - >::assume_init(res) + , Const>>::assume_init(res) } } diff --git a/src/base/edition.rs b/src/base/edition.rs index 81e10b48..f013ffd3 100644 --- a/src/base/edition.rs +++ b/src/base/edition.rs @@ -4,6 +4,7 @@ use std::cmp; use std::iter::ExactSizeIterator; #[cfg(any(feature = "std", feature = "alloc"))] use std::mem; +use std::mem::MaybeUninit; use std::ptr; use crate::base::allocator::{Allocator, Reallocator}; @@ -49,13 +50,10 @@ impl> Matrix { where I: IntoIterator, I::IntoIter: ExactSizeIterator + Clone, - 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 mut res = OMatrix::::new_uninitialized_generic(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. @@ -71,12 +69,12 @@ impl> Matrix { for (destination, source) in irows.clone().enumerate() { unsafe { *res.vget_unchecked_mut(destination) = - src.vget_unchecked(*source).inlined_clone() + MaybeUninit::new(src.vget_unchecked(*source).inlined_clone()); } } } - res + unsafe { res.assume_init() } } /// Creates a new matrix by extracting the given set of columns from `self`. @@ -90,15 +88,19 @@ impl> Matrix { { 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 mut res = Matrix::new_uninitialized_generic(nrows, Dynamic::new(icols.len())); for (destination, source) in icols.enumerate() { - res.column_mut(destination).copy_from(&self.column(*source)) + for (d, s) in res + .column_mut(destination) + .iter_mut() + .zip(self.column(*source).iter()) + { + *d = MaybeUninit::new(s.clone()); + } } - res + unsafe { res.assume_init() } } } @@ -190,7 +192,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: Clone, + { let (nrows, ncols) = self.shape(); let n = cmp::min(nrows, ncols); @@ -201,19 +206,25 @@ 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: Clone, + { 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() } + unsafe { *self.get_unchecked_mut((i, j)) = val.clone() } } } /// 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: Clone, + { 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() } + unsafe { *self.get_unchecked_mut((i, j)) = val.clone() } } } @@ -225,10 +236,13 @@ 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: Clone, + { for j in 0..self.ncols() { for i in (j + shift)..self.nrows() { - unsafe { *self.get_unchecked_mut((i, j)) = val.inlined_clone() } + unsafe { *self.get_unchecked_mut((i, j)) = val.clone() } } } } @@ -241,12 +255,15 @@ 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: Clone, + { for j in shift..self.ncols() { // TODO: is there a more efficient way to avoid the min ? // (necessary for rectangular matrices) for i in 0..cmp::min(j + 1 - shift, self.nrows()) { - unsafe { *self.get_unchecked_mut((i, j)) = val.inlined_clone() } + unsafe { *self.get_unchecked_mut((i, j)) = val.clone() } } } } @@ -921,9 +938,8 @@ impl OMatrix { where DefaultAllocator: Reallocator, { - let placeholder = unsafe { - crate::unimplemented_or_uninitialized_generic!(Dynamic::new(0), Dynamic::new(0)) - }; + let placeholder = + Matrix::new_uninitialized_generic(Dynamic::new(0), Dynamic::new(0)).assume_init(); let old = mem::replace(self, placeholder); let new = old.resize(new_nrows, new_ncols, val); let _ = mem::replace(self, new); @@ -946,9 +962,7 @@ where where DefaultAllocator: Reallocator, { - let placeholder = unsafe { - crate::unimplemented_or_uninitialized_generic!(Dynamic::new(0), self.data.shape().1) - }; + let placeholder = Matrix::from_fn_generic(Dynamic::new(0), self.data.shape().1, |_, _| val); let old = mem::replace(self, placeholder); let new = old.resize_vertically(new_nrows, val); let _ = mem::replace(self, new); @@ -971,9 +985,7 @@ where where DefaultAllocator: Reallocator, { - let placeholder = unsafe { - crate::unimplemented_or_uninitialized_generic!(self.data.shape().0, Dynamic::new(0)) - }; + let placeholder = Matrix::from_fn_generic(self.data.shape().0, Dynamic::new(0), |_, _| val); let old = mem::replace(self, placeholder); let new = old.resize_horizontally(new_ncols, val); let _ = mem::replace(self, new); diff --git a/src/base/matrix.rs b/src/base/matrix.rs index 7e8f79cc..51c8b945 100644 --- a/src/base/matrix.rs +++ b/src/base/matrix.rs @@ -29,7 +29,7 @@ use crate::base::storage::{ ContiguousStorage, ContiguousStorageMut, Owned, SameShapeStorage, Storage, StorageMut, }; use crate::base::{Const, DefaultAllocator, OMatrix, OVector, Scalar, Unit}; -use crate::{ArrayStorage, SMatrix, SimdComplexField}; +use crate::{ArrayStorage, MatrixSlice, MatrixSliceMut, SMatrix, SimdComplexField}; #[cfg(any(feature = "std", feature = "alloc"))] use crate::{DMatrix, DVector, Dynamic, VecStorage}; @@ -347,16 +347,13 @@ impl Matrix { } } -impl OMatrix, R, C> +impl OMatrix where DefaultAllocator: Allocator, { /// Allocates a matrix with the given number of rows and columns without initializing its content. - /// - /// Note: calling `Self::new_uninitialized_generic` is often **not** what you want to do. Consider - /// calling `Matrix::new_uninitialized_generic` instead. - pub fn new_uninitialized_generic(nrows: R, ncols: C) -> Self { - Self { + pub fn new_uninitialized_generic(nrows: R, ncols: C) -> OMatrix, R, C> { + OMatrix { data: >::allocate_uninitialized(nrows, ncols), _phantoms: PhantomData, } @@ -376,6 +373,24 @@ where } } +impl Matrix, R, C, S> { + /// Creates a full slice from `self` and assumes it to be initialized. + pub unsafe fn assume_init_ref(&self) -> MatrixSlice + where + S: Storage, R, C>, + { + self.full_slice().slice_assume_init() + } + + /// Creates a full mutable slice from `self` and assumes it to be initialized. + pub unsafe fn assume_init_mut(&mut self) -> MatrixSliceMut + where + S: StorageMut, R, C>, + { + self.full_slice_mut().slice_assume_init() + } +} + impl SMatrix { /// Creates a new statically-allocated matrix from the given [ArrayStorage]. /// @@ -428,6 +443,7 @@ impl> Matrix { /// Creates a new uninitialized matrix with the given uninitialized data pub unsafe fn from_uninitialized_data(data: MaybeUninit) -> MaybeUninit { + // BEEP BEEP this doesn't seem good let res: Matrix> = Matrix { data, _phantoms: PhantomData, @@ -493,6 +509,7 @@ impl> Matrix { /// let slice = mat.slice_with_steps((0, 0), (5, 3), (1, 2)); /// // The column strides is the number of steps (here 2) multiplied by the corresponding dimension. /// assert_eq!(mat.strides(), (1, 10)); + /// ``` #[inline] #[must_use] pub fn strides(&self) -> (usize, usize) { @@ -657,7 +674,7 @@ impl> Matrix { } } - unsafe { res.assume_init()} + unsafe { res.assume_init() } } /// Transposes `self` and store the result into `out`, which will become @@ -815,7 +832,7 @@ impl> Matrix { { let (nrows, ncols) = self.data.shape(); - let mut res = OMatrix::::new_uninitialized_generic(nrows, ncols); + let mut res = OMatrix::new_uninitialized_generic(nrows, ncols); assert_eq!( (nrows.value(), ncols.value()), @@ -1201,13 +1218,25 @@ impl> Matrix { } } - /// Fills this matrix with the content of another one. Both must have the same shape. + /// Fills this matrix with the content of another one via clones. Both must have the same shape. #[inline] pub fn copy_from(&mut self, other: &Matrix) where T: Clone, SB: Storage, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, + { + self.copy_from_fn(other, T::clone) + } + + /// Fills this matrix with the content of another one, after applying a function to + /// the references of the entries of the other matrix. Both must have the same shape. + #[inline] + pub fn copy_from_fn(&mut self, other: &Matrix, f: F) + where + SB: Storage, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, + F: FnMut(&U) -> T, { assert!( self.shape() == other.shape(), @@ -1217,19 +1246,68 @@ impl> Matrix { for j in 0..self.ncols() { for i in 0..self.nrows() { unsafe { - *self.get_unchecked_mut((i, j)) = other.get_unchecked((i, j)).clone(); + *self.get_unchecked_mut((i, j)) = f(other.get_unchecked((i, j))); } } } } - /// Fills this matrix with the content of the transpose another one. + /// Fills this matrix with the content of another one, after applying a function to + /// the entries of the other matrix. Both must have the same shape. + #[inline] + pub fn move_from(&mut self, other: Matrix) + where + SB: Storage, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, + { + self.move_from_fn(other, |e| e) + } + + /// Fills this matrix with the content of another one via moves. Both must have the same shape. + #[inline] + pub fn move_from_fn(&mut self, other: Matrix, f: F) + where + SB: Storage, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, + F: FnMut(U) -> T, + { + assert!( + self.shape() == other.shape(), + "Unable to move from a matrix with a different shape." + ); + + for j in 0..self.ncols() { + for i in 0..self.nrows() { + unsafe { + *self.get_unchecked_mut((i, j)) = f(*other.get_unchecked((i, j))); + } + } + } + } + + /// Fills this matrix with the content of the transpose another one via clones. #[inline] pub fn tr_copy_from(&mut self, other: &Matrix) where T: Clone, SB: Storage, ShapeConstraint: DimEq + SameNumberOfColumns, + { + self.tr_copy_from_fn(other, T::clone) + } + + /// Fills this matrix with the content of the transpose of another one, after applying + /// a function to the references of the entries of the other matrix. Both must have the + /// same shape. + #[inline] + pub fn tr_copy_from_fn( + &mut self, + other: &Matrix, + f: F, + ) where + SB: Storage, + ShapeConstraint: DimEq + SameNumberOfColumns, + F: FnMut(&U) -> T, { let (nrows, ncols) = self.shape(); assert!( @@ -1240,7 +1318,44 @@ impl> Matrix { for j in 0..ncols { for i in 0..nrows { unsafe { - *self.get_unchecked_mut((i, j)) = other.get_unchecked((j, i)).clone(); + *self.get_unchecked_mut((i, j)) = f(other.get_unchecked((j, i))); + } + } + } + } + + /// Fills this matrix with the content of the transpose another one via moves. + #[inline] + pub fn tr_move_from(&mut self, other: Matrix) + where + SB: Storage, + ShapeConstraint: DimEq + SameNumberOfColumns, + { + self.tr_move_from_fn(other, |e| e) + } + + /// Fills this matrix with the content of the transpose of another one, after applying + /// a function to the entries of the other matrix. Both must have the same shape. + #[inline] + pub fn tr_move_from_fn( + &mut self, + other: Matrix, + f: F, + ) where + SB: Storage, + ShapeConstraint: DimEq + SameNumberOfColumns, + F: FnMut(U) -> T, + { + let (nrows, ncols) = self.shape(); + assert!( + (ncols, nrows) == other.shape(), + "Unable to move from a matrix with incompatible shape." + ); + + for j in 0..ncols { + for i in 0..nrows { + unsafe { + *self.get_unchecked_mut((i, j)) = f(*other.get_unchecked((j, i))); } } } @@ -1316,11 +1431,9 @@ 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, R2, C2, SB>) + pub fn adjoint_to(&self, out: &mut Matrix, R2, C2, SB>) where - R2: Dim, - C2: Dim, - SB: StorageMut, + SB: StorageMut, R2, C2>, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { let (nrows, ncols) = self.shape(); @@ -1348,23 +1461,20 @@ impl> Matrix, { let (nrows, ncols) = self.data.shape(); + let mut res = OMatrix::new_uninitialized_generic(ncols, nrows); + self.adjoint_to(&mut res); - unsafe { - let mut res = OMatrix::new_uninitialized_generic(ncols, nrows); - self.adjoint_to(&mut res); - - res - } + unsafe { res.assume_init() } } /// Takes the conjugate and transposes `self` and store the result into `out`. #[deprecated(note = "Renamed `self.adjoint_to(out)`.")] #[inline] - pub fn conjugate_transpose_to(&self, out: &mut Matrix) - where - R2: Dim, - C2: Dim, - SB: StorageMut, + pub fn conjugate_transpose_to( + &self, + out: &mut Matrix, R2, C2, SB>, + ) where + SB: StorageMut, R2, C2>, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, { self.adjoint_to(out) @@ -1495,7 +1605,7 @@ impl> SquareMatrix { ); let dim = self.data.shape().0; - let mut res = OVector::::new_uninitialized_generic(dim, Const::<1>); + let mut res = OVector::new_uninitialized_generic(dim, Const::<1>); for i in 0..dim.value() { unsafe { @@ -1505,7 +1615,7 @@ impl> SquareMatrix { } // Safety: we have initialized all entries. - unsafe { Matrix::assume_init(res) } + unsafe { res.assume_init() } } /// Computes a trace of a square matrix, i.e., the sum of its diagonal elements. @@ -1630,13 +1740,12 @@ 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>) }; + let mut res = OVector::new_uninitialized_generic(hnrows, Const::<1>); res.generic_slice_mut((0, 0), self.data.shape()) - .copy_from(self); - res[(len, 0)] = element; + .copy_from_fn(self, |e| MaybeUninit::new(e.clone())); + res[(len, 0)] = MaybeUninit::new(element); - res + unsafe { res.assume_init() } } } @@ -1953,10 +2062,11 @@ impl(&self, b: &Matrix) -> MatrixCross + pub fn cross( + &self, + b: &Matrix, + ) -> MatrixCross where - R2: Dim, - C2: Dim, SB: Storage, DefaultAllocator: SameShapeAllocator, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, @@ -1974,8 +2084,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::new_uninitialized_generic(nrows, ncols); let ax = self.get_unchecked((0, 0)); let ay = self.get_unchecked((1, 0)); @@ -1985,22 +2094,27 @@ 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::new_uninitialized_generic(nrows, ncols); let ax = self.get_unchecked((0, 0)); let ay = self.get_unchecked((0, 1)); @@ -2010,14 +2124,20 @@ impl + SliceStorage<'a, MaybeUninit, R, C, RStride, CStride> +{ + pub unsafe fn assume_init(self) -> SliceStorage<'a, T, R, C, RStride, CStride> { + Self::from_raw_parts(self.ptr as *const T, self.shape, self.strides) + } +} + +impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> + SliceStorageMut<'a, MaybeUninit, R, C, RStride, CStride> +{ + pub unsafe fn assume_init(self) -> SliceStorageMut<'a, T, R, C, RStride, CStride> { + Self::from_raw_parts(self.ptr as *mut T, self.shape, self.strides) + } +} + unsafe impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> StorageMut for SliceStorageMut<'a, T, R, C, RStride, CStride> { @@ -242,10 +259,12 @@ unsafe impl<'a, T, R: Dim, CStride: Dim> ContiguousStorage for SliceStorage<'a, T, R, U1, U1, CStride> { } + unsafe impl<'a, T, R: Dim, CStride: Dim> ContiguousStorage for SliceStorageMut<'a, T, R, U1, U1, CStride> { } + unsafe impl<'a, T, R: Dim, CStride: Dim> ContiguousStorageMut for SliceStorageMut<'a, T, R, U1, U1, CStride> { @@ -255,10 +274,12 @@ unsafe impl<'a, T, R: DimName, C: Dim + IsNotStaticOne> ContiguousStorage { } + unsafe impl<'a, T, R: DimName, C: Dim + IsNotStaticOne> ContiguousStorage for SliceStorageMut<'a, T, R, C, U1, R> { } + unsafe impl<'a, T, R: DimName, C: Dim + IsNotStaticOne> ContiguousStorageMut for SliceStorageMut<'a, T, R, C, U1, R> { @@ -312,6 +333,7 @@ macro_rules! matrix_slice_impl( $fixed_slice_with_steps: ident, $generic_slice: ident, $generic_slice_with_steps: ident, + $full_slice: ident, $rows_range_pair: ident, $columns_range_pair: ident) => { /* @@ -370,7 +392,7 @@ macro_rules! matrix_slice_impl( pub fn $rows_generic($me: $Me, row_start: usize, nrows: RSlice) -> $MatrixSlice { - let my_shape = $me.data.shape(); + let my_shape = $me.data.shape(); $me.assert_slice_index((row_start, 0), (nrows.value(), my_shape.1.value()), (0, 0)); let shape = (nrows, my_shape.1); @@ -388,12 +410,12 @@ macro_rules! matrix_slice_impl( -> $MatrixSlice where RSlice: Dim { - let my_shape = $me.data.shape(); + let my_shape = $me.data.shape(); let my_strides = $me.data.strides(); $me.assert_slice_index((row_start, 0), (nrows.value(), my_shape.1.value()), (step, 0)); let strides = (Dynamic::new((step + 1) * my_strides.0.value()), my_strides.1); - let shape = (nrows, my_shape.1); + let shape = (nrows, my_shape.1); unsafe { let data = $SliceStorage::new_with_strides_unchecked($data, (row_start, 0), shape, strides); @@ -468,20 +490,19 @@ macro_rules! matrix_slice_impl( } } - /// Extracts from this matrix `ncols` columns skipping `step` columns. Both argument may /// or may not be values known at compile-time. #[inline] pub fn $columns_generic_with_step($me: $Me, first_col: usize, ncols: CSlice, step: usize) -> $MatrixSlice { - let my_shape = $me.data.shape(); + let my_shape = $me.data.shape(); let my_strides = $me.data.strides(); $me.assert_slice_index((0, first_col), (my_shape.0.value(), ncols.value()), (0, step)); let strides = (my_strides.0, Dynamic::new((step + 1) * my_strides.1.value())); - let shape = (my_shape.0, ncols); + let shape = (my_shape.0, ncols); unsafe { let data = $SliceStorage::new_with_strides_unchecked($data, (0, first_col), shape, strides); @@ -509,7 +530,6 @@ macro_rules! matrix_slice_impl( } } - /// Slices this matrix starting at its component `(start.0, start.1)` and with /// `(shape.0, shape.1)` components. Each row (resp. column) of the sliced matrix is /// separated by `steps.0` (resp. `steps.1`) ignored rows (resp. columns) of the @@ -550,11 +570,9 @@ macro_rules! matrix_slice_impl( /// Creates a slice that may or may not have a fixed size and stride. #[inline] - pub fn $generic_slice($me: $Me, start: (usize, usize), shape: (RSlice, CSlice)) + pub fn $generic_slice($me: $Me, start: (usize, usize), shape: (RSlice, CSlice)) -> $MatrixSlice - where RSlice: Dim, - CSlice: Dim { - + { $me.assert_slice_index(start, (shape.0.value(), shape.1.value()), (0, 0)); unsafe { @@ -585,6 +603,12 @@ macro_rules! matrix_slice_impl( } } + /// Returns a slice containing the entire matrix. + pub fn $full_slice($me: $Me) -> $MatrixSlice { + let (nrows, ncols) = $me.shape(); + $me.generic_slice((0, 0), (R::from_usize(nrows), C::from_usize(ncols))) + } + /* * * Splitting. @@ -697,6 +721,7 @@ impl> Matrix { fixed_slice_with_steps, generic_slice, generic_slice_with_steps, + full_slice, rows_range_pair, columns_range_pair); } @@ -727,10 +752,27 @@ impl> Matrix { fixed_slice_with_steps_mut, generic_slice_mut, generic_slice_with_steps_mut, + full_slice_mut, rows_range_pair_mut, columns_range_pair_mut); } +impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> + MatrixSlice<'a, MaybeUninit, R, C, RStride, CStride> +{ + pub unsafe fn slice_assume_init(self) -> MatrixSlice<'a, T, R, C, RStride, CStride> { + Matrix::from_data(self.data.assume_init()) + } +} + +impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> + MatrixSliceMut<'a, MaybeUninit, R, C, RStride, CStride> +{ + pub unsafe fn slice_assume_init(self) -> MatrixSliceMut<'a, T, R, C, RStride, CStride> { + Matrix::from_data(self.data.assume_init()) + } +} + /// A range with a size that may be known at compile-time. /// /// This may be: diff --git a/src/base/ops.rs b/src/base/ops.rs index 8da0249f..44b1c7c5 100644 --- a/src/base/ops.rs +++ b/src/base/ops.rs @@ -24,7 +24,7 @@ use crate::SimdComplexField; * Indexing. * */ -impl> Index for Matrix { +impl> Index for Matrix { type Output = T; #[inline] @@ -36,7 +36,6 @@ impl> Index for Matrix Index<(usize, usize)> for Matrix where - T: Scalar, S: Storage, { type Output = T; @@ -54,7 +53,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); @@ -64,7 +63,6 @@ impl> IndexMut for Matr impl IndexMut<(usize, usize)> for Matrix where - T: Scalar, S: StorageMut, { #[inline] @@ -139,15 +137,15 @@ macro_rules! componentwise_binop_impl( $TraitAssign: ident, $method_assign: ident, $method_assign_statically_unchecked: ident, $method_assign_statically_unchecked_rhs: ident; $method_to: ident, $method_to_statically_unchecked: ident) => { - impl> Matrix - where T: Scalar + $bound { - + where + T: Scalar + $bound + { /* * * Methods without dimension checking at compile-time. - * This is useful for code reuse because the sum representative system does not plays - * easily with static checks. + * This is useful for code reuse because the sum representative system does not play + * nicely with static checks. * */ #[inline] @@ -155,7 +153,7 @@ macro_rules! componentwise_binop_impl( &self, rhs: &Matrix, out: &mut Matrix, R3, C3, SC> ) where SB: Storage, - SC: StorageMut + StorageMut, R3, C3> + SC: StorageMut, R3, C3> { assert_eq!(self.shape(), rhs.shape(), "Matrix addition/subtraction dimensions mismatch."); assert_eq!(self.shape(), out.shape(), "Matrix addition/subtraction output dimensions mismatch."); @@ -184,13 +182,13 @@ macro_rules! componentwise_binop_impl( } } - #[inline] - fn $method_assign_statically_unchecked(&mut self, rhs: &Matrix) - where R2: Dim, - C2: Dim, - SA: StorageMut, - SB: Storage { + fn $method_assign_statically_unchecked( + &mut self, rhs: &Matrix + ) where + SA: StorageMut, + SB: Storage + { assert_eq!(self.shape(), rhs.shape(), "Matrix addition/subtraction dimensions mismatch."); // This is the most common case and should be deduced at compile-time. @@ -213,12 +211,12 @@ macro_rules! componentwise_binop_impl( } } - #[inline] - fn $method_assign_statically_unchecked_rhs(&self, rhs: &mut Matrix) - where R2: Dim, - C2: Dim, - SB: StorageMut { + fn $method_assign_statically_unchecked_rhs( + &self, rhs: &mut Matrix + ) where + SB: StorageMut + { assert_eq!(self.shape(), rhs.shape(), "Matrix addition/subtraction dimensions mismatch."); // This is the most common case and should be deduced at compile-time. @@ -253,14 +251,19 @@ macro_rules! componentwise_binop_impl( */ /// Equivalent to `self + rhs` but stores the result into `out` to avoid allocations. #[inline] - pub fn $method_to(&self, - rhs: &Matrix, - out: &mut Matrix) - where SB: Storage, - SC: StorageMut, - ShapeConstraint: SameNumberOfRows + SameNumberOfColumns + - SameNumberOfRows + SameNumberOfColumns { + pub fn $method_to( + &self, + rhs: &Matrix, + out: &mut Matrix, R3, C3, SC> + ) where + SB: Storage, + SC: StorageMut, R3, C3>, + ShapeConstraint: + SameNumberOfRows + + SameNumberOfColumns + + SameNumberOfRows + + SameNumberOfColumns + { self.$method_to_statically_unchecked(rhs, out) } } @@ -283,13 +286,14 @@ macro_rules! componentwise_binop_impl( } } - impl<'a, T, R1, C1, R2, C2, SA, SB> $Trait> for &'a Matrix - where R1: Dim, C1: Dim, R2: Dim, C2: Dim, - T: Scalar + $bound, - SA: Storage, - SB: Storage, - DefaultAllocator: SameShapeAllocator, - ShapeConstraint: SameNumberOfRows + SameNumberOfColumns { + impl<'a, T, R1: Dim, C1: Dim, R2: Dim, C2: Dim, SA, SB> $Trait> for &'a Matrix + where + T: Scalar + $bound, + SA: Storage, + SB: Storage, + DefaultAllocator: SameShapeAllocator, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns + { type Output = MatrixSum; #[inline] @@ -301,13 +305,14 @@ macro_rules! componentwise_binop_impl( } } - impl $Trait> for Matrix - where R1: Dim, C1: Dim, R2: Dim, C2: Dim, - T: Scalar + $bound, - SA: Storage, - SB: Storage, - DefaultAllocator: SameShapeAllocator, - ShapeConstraint: SameNumberOfRows + SameNumberOfColumns { + impl $Trait> for Matrix + where + T: Scalar + $bound, + SA: Storage, + SB: Storage, + DefaultAllocator: SameShapeAllocator, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns + { type Output = MatrixSum; #[inline] @@ -316,49 +321,48 @@ macro_rules! componentwise_binop_impl( } } - impl<'a, 'b, T, R1, C1, R2, C2, SA, SB> $Trait<&'b Matrix> for &'a Matrix - where R1: Dim, C1: Dim, R2: Dim, C2: Dim, - T: Scalar + $bound, - SA: Storage, - SB: Storage, - DefaultAllocator: SameShapeAllocator, - ShapeConstraint: SameNumberOfRows + SameNumberOfColumns { + impl<'a, 'b, T, R1: Dim, C1: Dim, R2: Dim, C2: Dim, SA, SB> $Trait<&'b Matrix> for &'a Matrix + where + T: Scalar + $bound, + SA: Storage, + SB: Storage, + DefaultAllocator: SameShapeAllocator, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns + { type Output = MatrixSum; #[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) - }; + let (nrows, ncols) = self.shape(); + let nrows: SameShapeR = Dim::from_usize(nrows); + let ncols: SameShapeC = Dim::from_usize(ncols); + let mut res = Matrix::new_uninitialized_generic(nrows, ncols); self.$method_to_statically_unchecked(rhs, &mut res); - res + unsafe { res.assume_init() } } } - impl<'b, T, R1, C1, R2, C2, SA, SB> $TraitAssign<&'b Matrix> for Matrix - where R1: Dim, C1: Dim, R2: Dim, C2: Dim, - T: Scalar + $bound, - SA: StorageMut, - SB: Storage, - ShapeConstraint: SameNumberOfRows + SameNumberOfColumns { - + impl<'b, T, R1: Dim, C1: Dim, R2: Dim, C2: Dim, SA, SB> $TraitAssign<&'b Matrix> for Matrix + where + T: Scalar + $bound, + SA: StorageMut, + SB: Storage, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns + { #[inline] fn $method_assign(&mut self, rhs: &'b Matrix) { self.$method_assign_statically_unchecked(rhs) } } - impl $TraitAssign> for Matrix - where R1: Dim, C1: Dim, R2: Dim, C2: Dim, - T: Scalar + $bound, - SA: StorageMut, - SB: Storage, - ShapeConstraint: SameNumberOfRows + SameNumberOfColumns { - + impl $TraitAssign> for Matrix + where + T: Scalar + $bound, + SA: StorageMut, + SB: Storage, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns + { #[inline] fn $method_assign(&mut self, rhs: Matrix) { self.$method_assign(&rhs) @@ -576,9 +580,9 @@ where #[inline] fn mul(self, rhs: &'b Matrix) -> Self::Output { - let mut res =Matrix::new_uninitialized_generic(self.data.shape().0, rhs.data.shape().1); - self.mul_to(rhs, &mut res); - unsafe{ res.assume_init()} + let mut res = Matrix::new_uninitialized_generic(self.data.shape().0, rhs.data.shape().1); + self.mul_to(rhs, &mut res); + unsafe { res.assume_init() } } } @@ -636,11 +640,8 @@ where // TODO: this is too restrictive: // − we can't use `a *= b` when `a` is a mutable slice. // − we can't use `a *= b` when C2 is not equal to C1. -impl MulAssign> for Matrix +impl MulAssign> for Matrix where - R1: Dim, - C1: Dim, - R2: Dim, T: Scalar + Zero + One + ClosedAdd + ClosedMul, SB: Storage, SA: ContiguousStorageMut + Clone, @@ -653,11 +654,8 @@ where } } -impl<'b, T, R1, C1, R2, SA, SB> MulAssign<&'b Matrix> for Matrix +impl<'b, T, R1: Dim, C1: Dim, R2: Dim, SA, SB> MulAssign<&'b Matrix> for Matrix where - R1: Dim, - C1: Dim, - R2: Dim, T: Scalar + Zero + One + ClosedAdd + ClosedMul, SB: Storage, SA: ContiguousStorageMut + Clone, @@ -697,7 +695,7 @@ where pub fn ad_mul(&self, rhs: &Matrix) -> OMatrix where T: SimdComplexField, - SB: Storage, R2, C2>, + SB: Storage, DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { @@ -746,7 +744,9 @@ 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)) = MaybeUninit::new(dot) ;} + unsafe { + *out.get_unchecked_mut((i, j)) = MaybeUninit::new(dot); + } } } } @@ -786,16 +786,16 @@ where #[inline] pub fn mul_to( &self, - rhs: &Matrix, R2, C2, SB>, - out: &mut Matrix, + rhs: &Matrix, + out: &mut Matrix, R3, C3, SC>, ) where SB: Storage, - SC: StorageMut, + SC: StorageMut, R3, C3>, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns + AreMultipliable, { - out.gemm(T::one(), self, rhs, T::zero()); + out.gemm_z(T::one(), self, rhs); } /// The kronecker product of two matrices (aka. tensor product of the corresponding linear diff --git a/src/base/statistics.rs b/src/base/statistics.rs index 59d78482..23ab524e 100644 --- a/src/base/statistics.rs +++ b/src/base/statistics.rs @@ -1,3 +1,5 @@ +use std::mem::MaybeUninit; + use crate::allocator::Allocator; use crate::storage::Storage; use crate::{Const, DefaultAllocator, Dim, Matrix, OVector, RowOVector, Scalar, VectorSlice, U1}; @@ -18,13 +20,12 @@ impl> Matrix { DefaultAllocator: Allocator, { let ncols = self.data.shape().1; - let mut res: RowOVector = - unsafe { crate::unimplemented_or_uninitialized_generic!(Const::<1>, ncols) }; + let mut res = RowOVector::new_uninitialized_generic(Const::<1>, ncols); for i in 0..ncols.value() { // TODO: avoid bound checking of column. unsafe { - *res.get_unchecked_mut((0, i)) = f(self.column(i)); + *res.get_unchecked_mut((0, i)) =MaybeUninit::new( f(self.column(i))); } } @@ -45,17 +46,16 @@ impl> Matrix { DefaultAllocator: Allocator, { let ncols = self.data.shape().1; - let mut res: OVector = - unsafe { crate::unimplemented_or_uninitialized_generic!(ncols, Const::<1>) }; + let mut res = Matrix::new_uninitialized_generic(ncols, Const::<1>); for i in 0..ncols.value() { // TODO: avoid bound checking of column. unsafe { - *res.vget_unchecked_mut(i) = f(self.column(i)); + *res.vget_unchecked_mut(i) = MaybeUninit::new(f(self.column(i))); } } - res + unsafe { res.assume_init() } } /// Returns a column vector resulting from the folding of `f` on each column of this matrix. diff --git a/src/linalg/permutation_sequence.rs b/src/linalg/permutation_sequence.rs index ea868b5a..a088c458 100644 --- a/src/linalg/permutation_sequence.rs +++ b/src/linalg/permutation_sequence.rs @@ -1,3 +1,5 @@ +use std::mem::MaybeUninit; + #[cfg(feature = "serde-serialize-no-std")] use serde::{Deserialize, Serialize}; @@ -8,7 +10,7 @@ use crate::allocator::Allocator; use crate::base::{DefaultAllocator, Matrix, OVector, Scalar}; #[cfg(any(feature = "std", feature = "alloc"))] use crate::dimension::Dynamic; -use crate::dimension::{Const, Dim, DimName}; +use crate::dimension::{ Dim, DimName}; use crate::storage::StorageMut; /// A sequence of row or column permutations. @@ -29,13 +31,13 @@ where DefaultAllocator: Allocator<(usize, usize), D>, { len: usize, - ipiv: OVector<(usize, usize), D>, + ipiv: OVector, D>, } impl Copy for PermutationSequence where DefaultAllocator: Allocator<(usize, usize), D>, - OVector<(usize, usize), D>: Copy, + OVector, D>: Copy, { } @@ -72,7 +74,7 @@ where unsafe { Self { len: 0, - ipiv: crate::unimplemented_or_uninitialized_generic!(dim, Const::<1>), + ipiv: OVector::new_uninitialized(dim), } } } @@ -97,7 +99,7 @@ where where S2: StorageMut, { - for i in self.ipiv.rows_range(..self.len).iter() { + for i in self.ipiv.rows_range(..self.len).iter().map(MaybeUninit::assume_init) { rhs.swap_rows(i.0, i.1) } }