From 9a528e23b9d14be126223532a069a621e8fe671b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Violeta=20Hern=C3=A1ndez?= Date: Sat, 17 Jul 2021 04:36:14 -0500 Subject: [PATCH] Almost! --- nalgebra-lapack/src/cholesky.rs | 6 +- src/base/blas.rs | 94 +++++++++++++++++++++++--- src/base/conversion.rs | 9 +-- src/base/default_allocator.rs | 22 +++--- src/base/edition.rs | 21 +++--- src/base/matrix.rs | 8 +-- src/base/statistics.rs | 3 +- src/geometry/dual_quaternion.rs | 3 +- src/geometry/orthographic.rs | 4 +- src/geometry/point.rs | 11 +-- src/geometry/point_conversion.rs | 3 +- src/geometry/transform.rs | 4 +- src/geometry/transform_ops.rs | 5 +- src/geometry/translation_conversion.rs | 9 ++- src/linalg/bidiagonal.rs | 14 ++-- src/linalg/cholesky.rs | 4 +- src/linalg/permutation_sequence.rs | 4 +- src/linalg/pow.rs | 8 ++- src/linalg/qr.rs | 57 +++++++++++++--- src/linalg/schur.rs | 92 ++++++++++++++++--------- src/linalg/svd.rs | 47 +++++++++++-- src/linalg/symmetric_eigen.rs | 42 ++++++++++-- src/linalg/symmetric_tridiagonal.rs | 57 ++++++++++++---- src/linalg/udu.rs | 41 +++++++++-- 24 files changed, 423 insertions(+), 145 deletions(-) diff --git a/nalgebra-lapack/src/cholesky.rs b/nalgebra-lapack/src/cholesky.rs index bc3515a5..929f2d40 100644 --- a/nalgebra-lapack/src/cholesky.rs +++ b/nalgebra-lapack/src/cholesky.rs @@ -24,17 +24,17 @@ use lapack; OMatrix: Deserialize<'de>")) )] #[derive(Clone, Debug)] -pub struct Cholesky +pub struct Cholesky where DefaultAllocator: Allocator, { l: OMatrix, } -impl Copy for Cholesky +impl Copy for Cholesky where DefaultAllocator: Allocator, - OMatrix: Copy, + Owned: Copy, { } diff --git a/src/base/blas.rs b/src/base/blas.rs index dec0af86..dd36ab37 100644 --- a/src/base/blas.rs +++ b/src/base/blas.rs @@ -329,22 +329,22 @@ where if !b.is_zero() { 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 { 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(); - } + } } } @@ -788,17 +788,89 @@ where for j in 1..ncols2 { let col2 = a.column(j); - let val = unsafe { x.vget_unchecked(j).inlined_clone() }; + let val = x.vget_unchecked(j).inlined_clone() ; init.axcpy(alpha.inlined_clone(), &col2, val, T::one()); } } } + + #[inline(always)] + fn xxgemv_z( + &mut self, + alpha: T, + a: &SquareMatrix, + x: &Vector, + dot: impl Fn( + &DVectorSlice, + &DVectorSlice, + ) -> T, + ) where + T: One, + SB: Storage, + SC: Storage, + ShapeConstraint: DimEq + AreMultipliable, + { + let dim1 = self.nrows(); + let dim2 = a.nrows(); + let dim3 = x.nrows(); + + assert!( + a.is_square(), + "Symmetric cgemv: the input matrix must be square." + ); + assert!( + dim2 == dim3 && dim1 == dim2, + "Symmetric cgemv: dimensions mismatch." + ); + + if dim2 == 0 { + return; + } + + // TODO: avoid bound checks. + let col2 = a.column(0); + let val = unsafe { x.vget_unchecked(0).inlined_clone() }; + self.axc(alpha.inlined_clone(), &col2, val); + + let mut res = unsafe { self.assume_init_mut() }; + res[0] += alpha.inlined_clone() * dot(&a.slice_range(1.., 0), &x.rows_range(1..)); + + for j in 1..dim2 { + let col2 = a.column(j); + let dot = dot(&col2.rows_range(j..), &x.rows_range(j..)); + + let val; + unsafe { + val = x.vget_unchecked(j).inlined_clone(); + *res.vget_unchecked_mut(j) += alpha.inlined_clone() * dot; + } + res.rows_range_mut(j + 1..).axpy( + alpha.inlined_clone() * val, + &col2.rows_range(j + 1..), + T::one(), + ); + } + } + + pub fn hegemv_z( + &mut self, + alpha: T, + a: &SquareMatrix, + x: &Vector, + ) where + T: SimdComplexField, + SB: Storage, + SC: Storage, + ShapeConstraint: DimEq + AreMultipliable, + { + self.xxgemv_z(alpha, a, x, |a, b| a.dotc(b)) + } } impl, R1, C1>> Matrix, R1, C1, S> where T: Scalar + Zero + One + ClosedAdd + ClosedMul, - // DefaultAllocator: Allocator, + // DefaultAllocator: Allocator, { /// Computes `alpha * a * b`, where `a` and `b` are matrices, and `alpha` is /// a scalar. @@ -850,7 +922,7 @@ where // 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(); + let (_, ncols3) = b.shape(); // Threshold determined empirically. const SMALL_DIM: usize = 5; @@ -1502,9 +1574,9 @@ where ShapeConstraint: DimEq + DimEq, DefaultAllocator: Allocator, { - let work = Matrix::new_uninitialized_generic(R3::from_usize(self.shape().0), Const::<1>); + let mut 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() }; + let mut work = unsafe { work.assume_init() }; self.ger(alpha.inlined_clone(), &work, &lhs.column(0), beta); @@ -1552,9 +1624,9 @@ where DefaultAllocator: Allocator, { // 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>); + let mut 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() }; + let mut work = unsafe { work.assume_init() }; self.column_mut(0) .gemv_tr(alpha.inlined_clone(), rhs, &work, beta.inlined_clone()); diff --git a/src/base/conversion.rs b/src/base/conversion.rs index f8e803fe..66ebe3bd 100644 --- a/src/base/conversion.rs +++ b/src/base/conversion.rs @@ -109,13 +109,14 @@ impl From<[T; D]> for SVector { } } -impl From> for [T; D] { +impl From> for [T; D] +where + T: Clone, +{ #[inline] fn from(vec: SVector) -> Self { // TODO: unfortunately, we must clone because we can move out of an array. - - // Counterpoint: this seems to work? - vec.data.0[0] + vec.data.0[0].clone() } } diff --git a/src/base/default_allocator.rs b/src/base/default_allocator.rs index 519f85f3..4551bcff 100644 --- a/src/base/default_allocator.rs +++ b/src/base/default_allocator.rs @@ -31,9 +31,9 @@ type DefaultUninitBuffer = * Allocator. * */ - /// A helper struct that controls how the storage for a matrix should be allocated. - /// - /// This struct is useless on its own. Instead, it's used in trait +/// A helper struct that controls how the storage for a matrix should be allocated. +/// +/// This struct is useless on its own. Instead, it's used in trait /// An allocator based on `GenericArray` and `VecStorage` for statically-sized and dynamically-sized /// matrices respectively. pub struct DefaultAllocator; @@ -72,7 +72,9 @@ impl Allocator, Const> for Def _: Const, _: Const, ) -> Owned, Const, Const> { - ArrayStorage([[MaybeUninit::uninit(); R]; C]) + // SAFETY: An uninitialized `[MaybeUninit<_>; LEN]` is valid. + let array = unsafe { MaybeUninit::uninit().assume_init() }; + ArrayStorage(array) } #[inline] @@ -126,9 +128,8 @@ impl Allocator for DefaultAllocator { let mut data = ManuallyDrop::new(uninit.data); // Safety: MaybeUninit has the same alignment and layout as T. - let new_data = unsafe { - Vec::from_raw_parts(data.as_mut_ptr() as *mut T, data.len(), data.capacity()) - }; + let new_data = + Vec::from_raw_parts(data.as_mut_ptr() as *mut T, data.len(), data.capacity()); VecStorage::new(uninit.nrows, uninit.ncols, new_data) } @@ -170,9 +171,8 @@ impl Allocator for DefaultAllocator { let mut data = ManuallyDrop::new(uninit.data); // Safety: MaybeUninit has the same alignment and layout as T. - let new_data = unsafe { - Vec::from_raw_parts(data.as_mut_ptr() as *mut T, data.len(), data.capacity()) - }; + let new_data = + Vec::from_raw_parts(data.as_mut_ptr() as *mut T, data.len(), data.capacity()); VecStorage::new(uninit.nrows, uninit.ncols, new_data) } @@ -184,7 +184,7 @@ impl Allocator for DefaultAllocator { * */ // Anything -> Static × Static -impl +impl Reallocator, Const> for DefaultAllocator where Self: Allocator, diff --git a/src/base/edition.rs b/src/base/edition.rs index 62977493..4e11bb26 100644 --- a/src/base/edition.rs +++ b/src/base/edition.rs @@ -178,7 +178,7 @@ impl> Matrix { /// Sets all the elements of this matrix to `f()`. #[inline] - pub fn fill_fn T>(&mut self, f: F) { + pub fn fill_fn T>(&mut self, mut f: F) { for e in self.iter_mut() { *e = f(); } @@ -942,8 +942,11 @@ impl OMatrix { where DefaultAllocator: Reallocator, { - let placeholder = - Matrix::new_uninitialized_generic(Dynamic::new(0), Dynamic::new(0)).assume_init(); + // BEEEP!!!! BEEEEEEEP!!! + + let placeholder = unsafe { + 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); @@ -966,7 +969,8 @@ where where DefaultAllocator: Reallocator, { - let placeholder = Matrix::from_fn_generic(Dynamic::new(0), self.data.shape().1, |_, _| val); + let placeholder = + Matrix::from_fn_generic(Dynamic::new(0), self.data.shape().1, |_, _| val.clone()); let old = mem::replace(self, placeholder); let new = old.resize_vertically(new_nrows, val); let _ = mem::replace(self, new); @@ -989,7 +993,8 @@ where where DefaultAllocator: Reallocator, { - let placeholder = Matrix::from_fn_generic(self.data.shape().0, Dynamic::new(0), |_, _| val); + let placeholder = + Matrix::from_fn_generic(self.data.shape().0, Dynamic::new(0), |_, _| val.clone()); let old = mem::replace(self, placeholder); let new = old.resize_horizontally(new_ncols, val); let _ = mem::replace(self, new); @@ -1059,11 +1064,7 @@ unsafe fn extend_rows(data: &mut [T], nrows: usize, ncols: usize, i: usize, n /// Extend the number of columns of the `Matrix` with elements from /// a given iterator. #[cfg(any(feature = "std", feature = "alloc"))] -impl Extend for Matrix -where - R: Dim, - S: Extend, -{ +impl> Extend for Matrix { /// Extend the number of columns of the `Matrix` with elements /// from the given iterator. /// diff --git a/src/base/matrix.rs b/src/base/matrix.rs index d13a467e..f973504b 100644 --- a/src/base/matrix.rs +++ b/src/base/matrix.rs @@ -1249,7 +1249,7 @@ impl> Matrix { /// 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) + pub fn copy_from_fn(&mut self, other: &Matrix,mut f: F) where SB: Storage, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, @@ -1282,7 +1282,7 @@ impl> Matrix { /// 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) + pub fn move_from_fn(&mut self, other: Matrix, mut f: F) where SB: Storage, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, @@ -1322,7 +1322,7 @@ impl> Matrix { pub fn tr_copy_from_fn( &mut self, other: &Matrix, - f: F, + mut f: F, ) where SB: Storage, ShapeConstraint: DimEq + SameNumberOfColumns, @@ -1359,7 +1359,7 @@ impl> Matrix { pub fn tr_move_from_fn( &mut self, other: Matrix, - f: F, + mut f: F, ) where SB: Storage, ShapeConstraint: DimEq + SameNumberOfColumns, diff --git a/src/base/statistics.rs b/src/base/statistics.rs index 2bb5ba7a..88f9236a 100644 --- a/src/base/statistics.rs +++ b/src/base/statistics.rs @@ -59,11 +59,12 @@ impl> Matrix { } /// Returns a column vector resulting from the folding of `f` on each column of this matrix. + // BEEEEP!!!! Pretty sure there's something fishy here. #[inline] #[must_use] pub fn compress_columns( &self, - init: OVector, + mut init: OVector, f: impl Fn(&mut OVector, VectorSlice), ) -> OVector where diff --git a/src/geometry/dual_quaternion.rs b/src/geometry/dual_quaternion.rs index 0fd10590..2c5968ef 100644 --- a/src/geometry/dual_quaternion.rs +++ b/src/geometry/dual_quaternion.rs @@ -278,7 +278,8 @@ impl<'a, T: Deserialize<'a>> Deserialize<'a> for DualQuaternion { } impl DualQuaternion { - fn to_vector(self) -> OVector { + // TODO: Cloning shouldn't be necessary. + fn to_vector(self) -> OVectorwhere T:Clone { (*self.as_ref()).into() } } diff --git a/src/geometry/orthographic.rs b/src/geometry/orthographic.rs index 98fd6b0d..974df3ff 100644 --- a/src/geometry/orthographic.rs +++ b/src/geometry/orthographic.rs @@ -28,7 +28,9 @@ impl Copy for Orthographic3 {} impl Clone for Orthographic3 { #[inline] fn clone(&self) -> Self { - Self::from_matrix_unchecked(self.matrix) + Self { + matrix: self.matrix.clone(), + } } } diff --git a/src/geometry/point.rs b/src/geometry/point.rs index a393bc2d..f65813e9 100644 --- a/src/geometry/point.rs +++ b/src/geometry/point.rs @@ -215,7 +215,7 @@ where let mut res = OVector::<_, DimNameSum>::new_uninitialized(); for i in 0..D::dim() { unsafe { - *res.get_unchecked(i) = MaybeUninit::new(self.coords[i].clone()); + *res.get_unchecked_mut(i) = MaybeUninit::new(self.coords[i].clone()); } } @@ -236,15 +236,16 @@ where // to avoid double-dropping. for i in 0..D::dim() { unsafe { - *res.get_unchecked(i) = MaybeUninit::new(self.coords[i]); + *res.get_unchecked_mut(i) = MaybeUninit::new(*self.coords.get_unchecked(i)); } } // Fix double drop - res[(D::dim(), 0)] = MaybeUninit::new(T::one()); - - unsafe { res.assume_init() } + unsafe { + *res.get_unchecked_mut(D::dim()) = MaybeUninit::new(T::one()); + res.assume_init() + } } /// Creates a new point with the given coordinates. diff --git a/src/geometry/point_conversion.rs b/src/geometry/point_conversion.rs index 022a7bd4..02ca1895 100644 --- a/src/geometry/point_conversion.rs +++ b/src/geometry/point_conversion.rs @@ -91,7 +91,8 @@ impl From<[T; D]> for Point { } } -impl From> for [T; D] { +impl From> for [T; D] where +T: Clone,{ #[inline] fn from(p: Point) -> Self { p.coords.into() diff --git a/src/geometry/transform.rs b/src/geometry/transform.rs index 1607a0b0..14bd43ae 100755 --- a/src/geometry/transform.rs +++ b/src/geometry/transform.rs @@ -550,8 +550,8 @@ where Const: DimNameAdd, C: SubTCategoryOf, DefaultAllocator: Allocator, U1>, DimNameSum, U1>> - + Allocator, U1>>, // + Allocator - // + Allocator + + Allocator, U1>>, + Owned, U1>, DimNameSum, U1>>: Clone, { /// Transform the given point by the inverse of this transformation. /// This may be cheaper than inverting the transformation and transforming diff --git a/src/geometry/transform_ops.rs b/src/geometry/transform_ops.rs index c4ec5cfc..8a21afd0 100644 --- a/src/geometry/transform_ops.rs +++ b/src/geometry/transform_ops.rs @@ -8,7 +8,7 @@ use simba::scalar::{ClosedAdd, ClosedMul, RealField, SubsetOf}; use crate::base::allocator::Allocator; use crate::base::dimension::{DimNameAdd, DimNameSum, U1}; -use crate::base::{Const, DefaultAllocator, OMatrix, SVector, Scalar}; +use crate::base::{Const, DefaultAllocator, OMatrix, SVector, Scalar};use crate::storage::Owned; use crate::geometry::{ Isometry, Point, Rotation, Similarity, SubTCategoryOf, SuperTCategoryOf, TAffine, TCategory, @@ -586,7 +586,8 @@ md_assign_impl_all!( const D; for CA, CB; where Const: DimNameAdd, CA: SuperTCategoryOf, CB: SubTCategoryOf, - DefaultAllocator: Allocator, U1>, DimNameSum, U1>>; + DefaultAllocator: Allocator, U1>, DimNameSum, U1>>, + Owned, U1>, DimNameSum, U1>>: Clone; self: Transform, rhs: Transform; [val] => #[allow(clippy::suspicious_op_assign_impl)] { *self *= rhs.inverse() }; [ref] => #[allow(clippy::suspicious_op_assign_impl)] { *self *= rhs.clone().inverse() }; diff --git a/src/geometry/translation_conversion.rs b/src/geometry/translation_conversion.rs index 7c75d379..bed39f7a 100644 --- a/src/geometry/translation_conversion.rs +++ b/src/geometry/translation_conversion.rs @@ -26,8 +26,8 @@ use crate::Point; */ impl SubsetOf> for Translation -where - T2: SupersetOf, +where + T2: SupersetOf, { #[inline] fn to_superset(&self) -> Translation { @@ -215,7 +215,10 @@ impl From> for Translation { } } -impl From> for [T; D] { +impl From> for [T; D] +where + T: Clone, +{ #[inline] fn from(t: Translation) -> Self { t.vector.into() diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs index 46bb9029..f25981a2 100644 --- a/src/linalg/bidiagonal.rs +++ b/src/linalg/bidiagonal.rs @@ -185,11 +185,13 @@ where ); } - Bidiagonal { - uv: matrix, - diagonal: diagonal.assume_init(), - off_diagonal: off_diagonal.assume_init(), - upper_diagonal, + unsafe { + Bidiagonal { + uv: matrix, + diagonal: diagonal.assume_init(), + off_diagonal: off_diagonal.assume_init(), + upper_diagonal, + } } } @@ -300,7 +302,7 @@ where let axis = self.uv.slice_range(i, i + shift..); let mut axis_packed = axis_packed.rows_range_mut(i + shift..); axis_packed.tr_copy_init_from(&axis); - let mut axis_packed = unsafe { axis_packed.slice_assume_init() }; + let axis_packed = unsafe { axis_packed.slice_assume_init() }; // TODO: sometimes, the axis might have a zero magnitude. let refl = Reflection::new(Unit::new_unchecked(axis_packed), T::zero()); diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 375ae521..afd90c0a 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -307,7 +307,7 @@ where ); chol.slice_range_mut(j + 1.., j).copy_init_from(&new_colj); - let chol = unsafe { chol.assume_init() }; + let mut chol = unsafe { chol.assume_init() }; // update the bottom right corner let mut bottom_right_corner = chol.slice_range_mut(j + 1.., j + 1..); @@ -348,7 +348,7 @@ where .copy_init_from(&self.chol.slice_range(j + 1.., ..j)); chol.slice_range_mut(j.., j..) .copy_init_from(&self.chol.slice_range(j + 1.., j + 1..)); - let chol = unsafe { chol.assume_init() }; + let mut chol = unsafe { chol.assume_init() }; // updates the bottom right corner let mut bottom_right_corner = chol.slice_range_mut(j.., j..); diff --git a/src/linalg/permutation_sequence.rs b/src/linalg/permutation_sequence.rs index e4594520..2cdfdd41 100644 --- a/src/linalg/permutation_sequence.rs +++ b/src/linalg/permutation_sequence.rs @@ -99,11 +99,11 @@ where /// Creates a new sequence of D identity permutations. #[inline] pub fn identity_generic(dim: D) -> Self { - unsafe { + Self { len: 0, ipiv: OVector::new_uninitialized_generic(dim, Const::<1>), - } + } } diff --git a/src/linalg/pow.rs b/src/linalg/pow.rs index df513643..68eb9682 100644 --- a/src/linalg/pow.rs +++ b/src/linalg/pow.rs @@ -40,18 +40,24 @@ where // We use the buffer to hold the result of multiplier ^ 2, thus avoiding // extra allocations. + let (nrows, ncols) = self.data.shape(); let mut multiplier = self.clone_owned(); - let mut buf = self.clone_owned(); + + // TODO: ACTUALLY MAKE BUF USEFUL! BEEEEEEEEP!! // Exponentiation by squares. loop { if e % two == one { + let mut buf = Matrix::new_uninitialized_generic(nrows, ncols); self.mul_to(&multiplier, &mut buf); + let buf = unsafe { buf.assume_init() }; self.copy_from(&buf); } e /= two; + let mut buf = Matrix::new_uninitialized_generic(nrows, ncols); multiplier.mul_to(&multiplier, &mut buf); + let buf = unsafe { buf.assume_init() }; multiplier.copy_from(&buf); if e == zero { diff --git a/src/linalg/qr.rs b/src/linalg/qr.rs index 4bdbb364..4b7d919c 100644 --- a/src/linalg/qr.rs +++ b/src/linalg/qr.rs @@ -1,3 +1,5 @@ +use std::fmt; + use num::Zero; #[cfg(feature = "serde-serialize-no-std")] use serde::{Deserialize, Serialize}; @@ -6,7 +8,7 @@ use crate::allocator::{Allocator, Reallocator}; use crate::base::{DefaultAllocator, Matrix, OMatrix, OVector, Unit}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::dimension::{Const, Dim, DimMin, DimMinimum}; -use crate::storage::{Storage, StorageMut}; +use crate::storage::{Owned, Storage, StorageMut}; use simba::scalar::ComplexField; use crate::geometry::Reflection; @@ -28,8 +30,8 @@ use crate::linalg::householder; OMatrix: Deserialize<'de>, OVector>: Deserialize<'de>")) )] -#[derive(Clone, Debug)] -pub struct QR, C: Dim> + +pub struct QR, C: Dim> where DefaultAllocator: Allocator + Allocator>, { @@ -37,14 +39,42 @@ where diag: OVector>, } -impl, C: Dim> Copy for QR +impl, C: Dim> Copy for QR where DefaultAllocator: Allocator + Allocator>, - OMatrix: Copy, - OVector>: Copy, + Owned: Copy, + Owned>: Copy, { } +impl, C: Dim> Clone for QR +where + DefaultAllocator: Allocator + Allocator>, + Owned: Clone, + Owned>: Clone, +{ + fn clone(&self) -> Self { + Self { + qr: self.qr.clone(), + diag: self.diag.clone(), + } + } +} + +impl, C: Dim> fmt::Debug for QR +where + DefaultAllocator: Allocator + Allocator>, + Owned: fmt::Debug, + Owned>: fmt::Debug, +{ + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + f.debug_struct("QR") + .field("qr", &self.qr) + .field("diag", &self.diag) + .finish() + } +} + impl, C: Dim> QR where DefaultAllocator: Allocator + Allocator + Allocator>, @@ -54,18 +84,23 @@ where let (nrows, ncols) = matrix.data.shape(); let min_nrows_ncols = nrows.min(ncols); - let mut diag = - unsafe { crate::unimplemented_or_uninitialized_generic!(min_nrows_ncols, Const::<1>) }; + let mut diag = Matrix::new_uninitialized_generic(min_nrows_ncols, Const::<1>); if min_nrows_ncols.value() == 0 { - return QR { qr: matrix, diag }; + return Self { + qr: matrix, + diag: unsafe { diag.assume_init() }, + }; } for i in 0..min_nrows_ncols.value() { - householder::clear_column_unchecked(&mut matrix, &mut diag[i], i, 0, None); + householder::clear_column_unchecked(&mut matrix, diag[i].as_mut_ptr(), i, 0, None); } - QR { qr: matrix, diag } + Self { + qr: matrix, + diag: unsafe { diag.assume_init() }, + } } /// Retrieves the upper trapezoidal submatrix `R` of this decomposition. diff --git a/src/linalg/schur.rs b/src/linalg/schur.rs index f359900d..f93aec1e 100644 --- a/src/linalg/schur.rs +++ b/src/linalg/schur.rs @@ -1,16 +1,18 @@ #![allow(clippy::suspicious_operation_groupings)] +use std::cmp; +use std::fmt; +use std::mem::MaybeUninit; + #[cfg(feature = "serde-serialize-no-std")] use serde::{Deserialize, Serialize}; use approx::AbsDiffEq; use num_complex::Complex as NumComplex; use simba::scalar::{ComplexField, RealField}; -use std::cmp; -use std::mem::MaybeUninit; use crate::allocator::Allocator; use crate::base::dimension::{Const, Dim, DimDiff, DimSub, Dynamic, U1, U2}; -use crate::base::storage::Storage; +use crate::base::storage::{Owned, Storage}; use crate::base::{DefaultAllocator, OMatrix, OVector, SquareMatrix, Unit, Vector2, Vector3}; use crate::geometry::Reflection; @@ -32,8 +34,7 @@ use crate::linalg::Hessenberg; serde(bound(deserialize = "DefaultAllocator: Allocator, OMatrix: Deserialize<'de>")) )] -#[derive(Clone, Debug)] -pub struct Schur +pub struct Schur where DefaultAllocator: Allocator, { @@ -41,13 +42,39 @@ where t: OMatrix, } -impl Copy for Schur +impl Copy for Schur where DefaultAllocator: Allocator, - OMatrix: Copy, + Owned: Copy, { } +impl Clone for Schur +where + DefaultAllocator: Allocator, + Owned: Clone, +{ + fn clone(&self) -> Self { + Self { + q: self.q.clone(), + t: self.t.clone(), + } + } +} + +impl fmt::Debug for Schur +where + DefaultAllocator: Allocator, + Owned: fmt::Debug, +{ + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + f.debug_struct("Schur") + .field("q", &self.q) + .field("t", &self.t) + .finish() + } +} + impl Schur where D: DimSub, // For Hessenberg. @@ -73,8 +100,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 = OVector::new_uninitialized_generic(m.data.shape().0, Const::<1>); Self::do_decompose(m, &mut work, eps, max_niter, true) .map(|(q, t)| Schur { q: q.unwrap(), t }) @@ -82,7 +108,7 @@ where fn do_decompose( mut m: OMatrix, - work: &mut OVector, + work: &mut OVector, D>, eps: T::RealField, max_niter: usize, compute_q: bool, @@ -271,7 +297,9 @@ where } /// Computes the eigenvalues of the decomposed matrix. - fn do_eigenvalues(t: &OMatrix, out: &mut OVector) -> bool { + fn do_eigenvalues(t: &OMatrix, out: &mut OVector, D>) -> bool { + // TODO: check dropping stuff. + let dim = t.nrows(); let mut m = 0; @@ -279,7 +307,7 @@ where let n = m + 1; if t[(n, m)].is_zero() { - out[m] = t[(m, m)]; + out[m] = MaybeUninit::new(t[(m, m)]); m += 1; } else { // Complex eigenvalue. @@ -288,7 +316,7 @@ where } if m == dim - 1 { - out[m] = t[(m, m)]; + out[m] = MaybeUninit::new(t[(m, m)]); } true @@ -297,11 +325,13 @@ where /// Computes the complex eigenvalues of the decomposed matrix. fn do_complex_eigenvalues( t: &OMatrix, - out: &mut OVector, D>, + out: &mut OVector>, D>, ) where T: RealField, DefaultAllocator: Allocator, D>, { + // TODO: check for dropping behavior. + let dim = t.nrows(); let mut m = 0; @@ -309,7 +339,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. @@ -391,11 +421,9 @@ 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 = OVector::new_uninitialized_generic(self.t.data.shape().0, Const::<1>); if Self::do_eigenvalues(&self.t, &mut out) { - Some(out) + Some(unsafe { out.assume_init() }) } else { None } @@ -408,11 +436,9 @@ 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 = OVector::new_uninitialized_generic(self.t.data.shape().0, Const::<1>); Self::do_complex_eigenvalues(&self.t, &mut out); - out + unsafe { out.assume_init() } } } @@ -517,14 +543,14 @@ where /// Computes the eigenvalues of this matrix. #[must_use] pub fn eigenvalues(&self) -> Option> { + // TODO: check drop stuff. + assert!( self.is_square(), "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 = OVector::new_uninitialized_generic(self.data.shape().0, Const::<1>); // Special case for 2x2 matrices. if self.nrows() == 2 { @@ -533,9 +559,9 @@ where let me = self.fixed_slice::<2, 2>(0, 0); return match compute_2x2_eigvals(&me) { Some((a, b)) => { - work[0] = a; - work[1] = b; - Some(work) + work[0] = MaybeUninit::new(a); + work[1] = MaybeUninit::new(b); + Some(unsafe { work.assume_init() }) } None => None, }; @@ -551,7 +577,7 @@ where ) .unwrap(); if Schur::do_eigenvalues(&schur.1, &mut work) { - Some(work) + Some(unsafe { work.assume_init() }) } else { None } @@ -566,7 +592,7 @@ where DefaultAllocator: Allocator, D>, { let dim = self.data.shape().0; - let mut work = unsafe { crate::unimplemented_or_uninitialized_generic!(dim, Const::<1>) }; + let mut work = OVector::new_uninitialized_generic(dim, Const::<1>); let schur = Schur::do_decompose( self.clone_owned(), @@ -576,8 +602,8 @@ where false, ) .unwrap(); - let mut eig = unsafe { crate::unimplemented_or_uninitialized_generic!(dim, Const::<1>) }; + let mut eig = OVector::new_uninitialized_generic(dim, Const::<1>); Schur::do_complex_eigenvalues(&schur.1, &mut eig); - eig + unsafe { eig.assume_init() } } } diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 241f00ce..c8cf5501 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -1,3 +1,5 @@ +use std::fmt; + #[cfg(feature = "serde-serialize-no-std")] use serde::{Deserialize, Serialize}; @@ -8,7 +10,7 @@ use crate::allocator::Allocator; use crate::base::{DefaultAllocator, Matrix, Matrix2x3, OMatrix, OVector, Vector2}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::dimension::{Dim, DimDiff, DimMin, DimMinimum, DimSub, U1}; -use crate::storage::Storage; +use crate::storage::{Owned, Storage}; use simba::scalar::{ComplexField, RealField}; use crate::linalg::givens::GivensRotation; @@ -39,7 +41,6 @@ use crate::linalg::Bidiagonal; OVector>: Deserialize<'de>" )) )] -#[derive(Clone, Debug)] pub struct SVD, C: Dim> where DefaultAllocator: Allocator, C> @@ -59,12 +60,48 @@ where DefaultAllocator: Allocator, C> + Allocator> + Allocator>, - OMatrix>: Copy, - OMatrix, C>: Copy, - OVector>: Copy, + Owned>: Copy, + Owned, C>: Copy, + Owned>: Copy, { } +impl, C: Dim> Clone for SVD +where + DefaultAllocator: Allocator, C> + + Allocator> + + Allocator>, + Owned>: Clone, + Owned, C>: Clone, + Owned>: Clone, +{ + fn clone(&self) -> Self { + Self { + u: self.u.clone(), + v_t: self.v_t.clone(), + singular_values: self.singular_values.clone(), + } + } +} + +impl, C: Dim> fmt::Debug for SVD +where + DefaultAllocator: Allocator, C> + + Allocator> + + Allocator>, + Owned>: fmt::Debug, + Owned, C>: fmt::Debug, + Owned>: fmt::Debug, +{ + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + f.debug_struct("SVD") + .field("u", &self.u) + .field("v_t", &self.v_t) + .field("singular_values", &self.singular_values) + .finish() + } +} + impl, C: Dim> SVD where DimMinimum: DimSub, // for Bidiagonal. diff --git a/src/linalg/symmetric_eigen.rs b/src/linalg/symmetric_eigen.rs index 5ac6d5da..ad4d6be4 100644 --- a/src/linalg/symmetric_eigen.rs +++ b/src/linalg/symmetric_eigen.rs @@ -1,3 +1,5 @@ +use std::fmt; + #[cfg(feature = "serde-serialize-no-std")] use serde::{Deserialize, Serialize}; @@ -7,7 +9,7 @@ use num::Zero; use crate::allocator::Allocator; use crate::base::{DefaultAllocator, Matrix2, OMatrix, OVector, SquareMatrix, Vector2}; use crate::dimension::{Dim, DimDiff, DimSub, U1}; -use crate::storage::Storage; +use crate::storage::{Owned, Storage}; use simba::scalar::ComplexField; use crate::linalg::givens::GivensRotation; @@ -29,7 +31,6 @@ use crate::linalg::SymmetricTridiagonal; OVector: Deserialize<'de>, OMatrix: Deserialize<'de>")) )] -#[derive(Clone, Debug)] pub struct SymmetricEigen where DefaultAllocator: Allocator + Allocator, @@ -44,11 +45,39 @@ where impl Copy for SymmetricEigen where DefaultAllocator: Allocator + Allocator, - OMatrix: Copy, - OVector: Copy, + Owned: Copy, + Owned: Copy, { } +impl Clone for SymmetricEigen +where + DefaultAllocator: Allocator + Allocator, + Owned: Clone, + Owned: Clone, +{ + fn clone(&self) -> Self { + Self { + eigenvectors: self.eigenvectors.clone(), + eigenvalues: self.eigenvalues.clone(), + } + } +} + +impl fmt::Debug for SymmetricEigen +where + DefaultAllocator: Allocator + Allocator, + Owned: fmt::Debug, + Owned: fmt::Debug, +{ + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + f.debug_struct("SymmetricEigen") + .field("eigenvectors", &self.eigenvectors) + .field("eigenvalues", &self.eigenvalues) + .finish() + } +} + impl SymmetricEigen where DefaultAllocator: Allocator + Allocator, @@ -270,7 +299,10 @@ where /// /// This is useful if some of the eigenvalues have been manually modified. #[must_use] - pub fn recompose(&self) -> OMatrix { + pub fn recompose(&self) -> OMatrix + where + Owned: Clone, + { let mut u_t = self.eigenvectors.clone(); for i in 0..self.eigenvalues.len() { let val = self.eigenvalues[i]; diff --git a/src/linalg/symmetric_tridiagonal.rs b/src/linalg/symmetric_tridiagonal.rs index c7e87ba8..cff9dc11 100644 --- a/src/linalg/symmetric_tridiagonal.rs +++ b/src/linalg/symmetric_tridiagonal.rs @@ -1,10 +1,13 @@ +use std::fmt; +use std::mem::MaybeUninit; + #[cfg(feature = "serde-serialize-no-std")] 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 crate::storage::{Owned, Storage}; use simba::scalar::ComplexField; use crate::linalg::householder; @@ -25,8 +28,7 @@ use crate::linalg::householder; OMatrix: Deserialize<'de>, OVector>: Deserialize<'de>")) )] -#[derive(Clone, Debug)] -pub struct SymmetricTridiagonal> +pub struct SymmetricTridiagonal> where DefaultAllocator: Allocator + Allocator>, { @@ -34,14 +36,42 @@ where off_diagonal: OVector>, } -impl> Copy for SymmetricTridiagonal +impl> Copy for SymmetricTridiagonal where DefaultAllocator: Allocator + Allocator>, - OMatrix: Copy, - OVector>: Copy, + Owned: Copy, + Owned>: Copy, { } +impl> Clone for SymmetricTridiagonal +where + DefaultAllocator: Allocator + Allocator>, + Owned: Clone, + Owned>: Clone, +{ + fn clone(&self) -> Self { + Self { + tri: self.tri.clone(), + off_diagonal: self.off_diagonal.clone(), + } + } +} + +impl> fmt::Debug for SymmetricTridiagonal +where + DefaultAllocator: Allocator + Allocator>, + Owned: fmt::Debug, + Owned>: fmt::Debug, +{ + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + f.debug_struct("SymmetricTridiagonal") + .field("tri", &self.tri) + .field("off_diagonal", &self.off_diagonal) + .finish() + } +} + impl> SymmetricTridiagonal where DefaultAllocator: Allocator + Allocator>, @@ -61,24 +91,21 @@ 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 = OVector::new_uninitialized_generic(dim.sub(Const::<1>), Const::<1>); + let mut p = OVector::new_uninitialized_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..); - p.hegemv(crate::convert(2.0), &m, &axis, T::zero()); + p.hegemv_z(crate::convert(2.0), &m, &axis); + let p = unsafe { p.slice_assume_init() }; let dot = axis.dotc(&p); m.hegerc(-T::one(), &p, &axis, T::one()); @@ -89,7 +116,7 @@ where Self { tri: m, - off_diagonal, + off_diagonal: unsafe { off_diagonal.assume_init() }, } } diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs index 7b4a9cc9..8e1b068f 100644 --- a/src/linalg/udu.rs +++ b/src/linalg/udu.rs @@ -1,10 +1,12 @@ +use std::fmt; + #[cfg(feature = "serde-serialize-no-std")] use serde::{Deserialize, Serialize}; use crate::allocator::Allocator; use crate::base::{Const, DefaultAllocator, OMatrix, OVector}; use crate::dimension::Dim; -use crate::storage::Storage; +use crate::storage::{Owned, Storage}; use simba::scalar::RealField; /// UDU factorization. @@ -19,8 +21,7 @@ use simba::scalar::RealField; deserialize = "OVector: Deserialize<'de>, OMatrix: Deserialize<'de>" )) )] -#[derive(Clone, Debug)] -pub struct UDU +pub struct UDU where DefaultAllocator: Allocator + Allocator, { @@ -30,14 +31,42 @@ where pub d: OVector, } -impl Copy for UDU +impl Copy for UDU where DefaultAllocator: Allocator + Allocator, - OVector: Copy, - OMatrix: Copy, + Owned: Copy, + Owned: Copy, { } +impl Clone for UDU +where + DefaultAllocator: Allocator + Allocator, + Owned: Clone, + Owned: Clone, +{ + fn clone(&self) -> Self { + Self { + u: self.u.clone(), + d: self.d.clone(), + } + } +} + +impl fmt::Debug for UDU +where + DefaultAllocator: Allocator + Allocator, + Owned: fmt::Debug, + Owned: fmt::Debug, +{ + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + f.debug_struct("UDU") + .field("u", &self.u) + .field("d", &self.d) + .finish() + } +} + impl UDU where DefaultAllocator: Allocator + Allocator,