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<T, D, D>: Deserialize<'de>"))
 )]
 #[derive(Clone, Debug)]
-pub struct Cholesky<T: Scalar, D: Dim>
+pub struct Cholesky<T, D: Dim>
 where
     DefaultAllocator: Allocator<T, D, D>,
 {
     l: OMatrix<T, D, D>,
 }
 
-impl<T: Scalar + Copy, D: Dim> Copy for Cholesky<T, D>
+impl<T: Copy, D: Dim> Copy for Cholesky<T, D>
 where
     DefaultAllocator: Allocator<T, D, D>,
-    OMatrix<T, D, D>: Copy,
+    Owned<T, D, D>: 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<D2: Dim, D3: Dim, SB, SC>(
+        &mut self,
+        alpha: T,
+        a: &SquareMatrix<T, D2, SB>,
+        x: &Vector<T, D3, SC>,
+        dot: impl Fn(
+            &DVectorSlice<T, SB::RStride, SB::CStride>,
+            &DVectorSlice<T, SC::RStride, SC::CStride>,
+        ) -> T,
+    ) where
+        T: One,
+        SB: Storage<T, D2, D2>,
+        SC: Storage<T, D3>,
+        ShapeConstraint: DimEq<D, D2> + AreMultipliable<D2, D2, D3, U1>,
+    {
+        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<D2: Dim, D3: Dim, SB, SC>(
+        &mut self,
+        alpha: T,
+        a: &SquareMatrix<T, D2, SB>,
+        x: &Vector<T, D3, SC>,
+    ) where
+        T: SimdComplexField,
+        SB: Storage<T, D2, D2>,
+        SC: Storage<T, D3>,
+        ShapeConstraint: DimEq<D, D2> + AreMultipliable<D2, D2, D3, U1>,
+    {
+        self.xxgemv_z(alpha, a, x, |a, b| a.dotc(b))
+    }
 }
 
 impl<T, R1: Dim, C1: Dim, S: StorageMut<MaybeUninit<T>, R1, C1>> Matrix<MaybeUninit<T>, R1, C1, S>
 where
     T: Scalar + Zero + One + ClosedAdd + ClosedMul,
-   // DefaultAllocator: Allocator<T, R1, C1>,
+    // DefaultAllocator: Allocator<T, R1, C1>,
 {
     /// 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<D1, R3> + DimEq<C3, D4>,
         DefaultAllocator: Allocator<T, R3>,
     {
-        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<T, D3>,
     {
         // 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::<D3, D3, R4, S3, _>(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<T, const D: usize> From<[T; D]> for SVector<T, D> {
     }
 }
 
-impl<T, const D: usize> From<SVector<T, D>> for [T; D] {
+impl<T, const D: usize> From<SVector<T, D>> for [T; D]
+where
+    T: Clone,
+{
     #[inline]
     fn from(vec: SVector<T, D>) -> 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<T, R, C> =
  * 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<T, const R: usize, const C: usize> Allocator<T, Const<R>, Const<C>> for Def
         _: Const<R>,
         _: Const<C>,
     ) -> Owned<MaybeUninit<T>, Const<R>, Const<C>> {
-        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<T, C: Dim> Allocator<T, Dynamic, C> for DefaultAllocator {
         let mut data = ManuallyDrop::new(uninit.data);
 
         // Safety: MaybeUninit<T> 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<T, R: DimName> Allocator<T, R, Dynamic> for DefaultAllocator {
         let mut data = ManuallyDrop::new(uninit.data);
 
         // Safety: MaybeUninit<T> 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<T, R: DimName> Allocator<T, R, Dynamic> for DefaultAllocator {
  *
  */
 // Anything -> Static × Static
-impl<T, RFrom:Dim, CFrom:Dim, const RTO: usize, const CTO: usize>
+impl<T, RFrom: Dim, CFrom: Dim, const RTO: usize, const CTO: usize>
     Reallocator<T, RFrom, CFrom, Const<RTO>, Const<CTO>> for DefaultAllocator
 where
     Self: Allocator<T, RFrom, CFrom>,
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<T, R: Dim, C: Dim, S: StorageMut<T, R, C>> Matrix<T, R, C, S> {
 
     /// Sets all the elements of this matrix to `f()`.
     #[inline]
-    pub fn fill_fn<F: FnMut() -> T>(&mut self, f: F) {
+    pub fn fill_fn<F: FnMut() -> T>(&mut self, mut f: F) {
         for e in self.iter_mut() {
             *e = f();
         }
@@ -942,8 +942,11 @@ impl<T: Clone> OMatrix<T, Dynamic, Dynamic> {
     where
         DefaultAllocator: Reallocator<T, Dynamic, Dynamic, Dynamic, Dynamic>,
     {
-        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<T, Dynamic, C, Dynamic, C>,
     {
-        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<T, R, Dynamic, R, Dynamic>,
     {
-        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<T>(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<T, R, S> Extend<T> for Matrix<T, R, Dynamic, S>
-where
-    R: Dim,
-    S: Extend<T>,
-{
+impl<T, R: Dim, S: Extend<T>> Extend<T> for Matrix<T, R, Dynamic, S> {
     /// 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<T, R: Dim, C: Dim, S: StorageMut<T, R, C>> Matrix<T, R, C, S> {
     /// 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<U, R2: Dim, C2: Dim, SB, F>(&mut self, other: &Matrix<U, R2, C2, SB>, f: F)
+    pub fn copy_from_fn<U, R2: Dim, C2: Dim, SB, F>(&mut self, other: &Matrix<U, R2, C2, SB>,mut f: F)
     where
         SB: Storage<U, R2, C2>,
         ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
@@ -1282,7 +1282,7 @@ impl<T, R: Dim, C: Dim, S: StorageMut<T, R, C>> Matrix<T, R, C, S> {
 
     /// Fills this matrix with the content of another one via moves. Both must have the same shape.
     #[inline]
-    pub fn move_from_fn<U, R2: Dim, C2: Dim, SB, F>(&mut self, other: Matrix<U, R2, C2, SB>, f: F)
+    pub fn move_from_fn<U, R2: Dim, C2: Dim, SB, F>(&mut self, other: Matrix<U, R2, C2, SB>, mut f: F)
     where
         SB: Storage<U, R2, C2>,
         ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
@@ -1322,7 +1322,7 @@ impl<T, R: Dim, C: Dim, S: StorageMut<T, R, C>> Matrix<T, R, C, S> {
     pub fn tr_copy_from_fn<U, R2: Dim, C2: Dim, SB, F>(
         &mut self,
         other: &Matrix<U, R2, C2, SB>,
-        f: F,
+   mut     f: F,
     ) where
         SB: Storage<U, R2, C2>,
         ShapeConstraint: DimEq<R, C2> + SameNumberOfColumns<C, R2>,
@@ -1359,7 +1359,7 @@ impl<T, R: Dim, C: Dim, S: StorageMut<T, R, C>> Matrix<T, R, C, S> {
     pub fn tr_move_from_fn<U, R2: Dim, C2: Dim, SB, F>(
         &mut self,
         other: Matrix<U, R2, C2, SB>,
-        f: F,
+      mut  f: F,
     ) where
         SB: Storage<U, R2, C2>,
         ShapeConstraint: DimEq<R, C2> + SameNumberOfColumns<C, R2>,
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<T, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     }
 
     /// 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<T, R>,
+        mut init: OVector<T, R>,
         f: impl Fn(&mut OVector<T, R>, VectorSlice<T, R, S::RStride, S::CStride>),
     ) -> OVector<T, R>
     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<T> {
 }
 
 impl<T> DualQuaternion<T> {
-    fn to_vector(self) -> OVector<T, U8> {
+    // TODO: Cloning shouldn't be necessary.
+    fn to_vector(self) -> OVector<T, U8>where 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<T: Copy> Copy for Orthographic3<T> {}
 impl<T: Clone> Clone for Orthographic3<T> {
     #[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<D, U1>>::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<T, const D: usize> From<[T; D]> for Point<T, D> {
     }
 }
 
-impl<T, const D: usize> From<Point<T, D>> for [T; D] {
+impl<T, const D: usize> From<Point<T, D>> for [T; D] where
+T: Clone,{
     #[inline]
     fn from(p: Point<T, D>) -> 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<D>: DimNameAdd<U1>,
     C: SubTCategoryOf<TProjective>,
     DefaultAllocator: Allocator<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>
-        + Allocator<T, DimNameSum<Const<D>, U1>>, // + Allocator<T, D, D>
-                                                  // + Allocator<T, D>
+        + Allocator<T, DimNameSum<Const<D>, U1>>,
+    Owned<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, 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<D>: DimNameAdd<U1>, CA: SuperTCategoryOf<CB>, CB: SubTCategoryOf<TProjective>,
-          DefaultAllocator: Allocator<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>;
+          DefaultAllocator: Allocator<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>,
+          Owned<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>: Clone;
     self: Transform<T, CA, D>, rhs: Transform<T, CB, D>;
     [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<T1, T2, const D: usize> SubsetOf<Translation<T2, D>> for Translation<T1, D>
-where 
-    T2:  SupersetOf<T1>,
+where
+    T2: SupersetOf<T1>,
 {
     #[inline]
     fn to_superset(&self) -> Translation<T2, D> {
@@ -215,7 +215,10 @@ impl<T, const D: usize> From<Point<T, D>> for Translation<T, D> {
     }
 }
 
-impl<T, const D: usize> From<Translation<T, D>> for [T; D] {
+impl<T, const D: usize> From<Translation<T, D>> for [T; D]
+where
+    T: Clone,
+{
     #[inline]
     fn from(t: Translation<T, D>) -> 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<T, R, C>: Deserialize<'de>,
          OVector<T, DimMinimum<R, C>>: Deserialize<'de>"))
 )]
-#[derive(Clone, Debug)]
-pub struct QR<T: ComplexField, R: DimMin<C>, C: Dim>
+
+pub struct QR<T, R: DimMin<C>, C: Dim>
 where
     DefaultAllocator: Allocator<T, R, C> + Allocator<T, DimMinimum<R, C>>,
 {
@@ -37,14 +39,42 @@ where
     diag: OVector<T, DimMinimum<R, C>>,
 }
 
-impl<T: ComplexField, R: DimMin<C>, C: Dim> Copy for QR<T, R, C>
+impl<T: Copy, R: DimMin<C>, C: Dim> Copy for QR<T, R, C>
 where
     DefaultAllocator: Allocator<T, R, C> + Allocator<T, DimMinimum<R, C>>,
-    OMatrix<T, R, C>: Copy,
-    OVector<T, DimMinimum<R, C>>: Copy,
+    Owned<T, R, C>: Copy,
+    Owned<T, DimMinimum<R, C>>: Copy,
 {
 }
 
+impl<T: Clone, R: DimMin<C>, C: Dim> Clone for QR<T, R, C>
+where
+    DefaultAllocator: Allocator<T, R, C> + Allocator<T, DimMinimum<R, C>>,
+    Owned<T, R, C>: Clone,
+    Owned<T, DimMinimum<R, C>>: Clone,
+{
+    fn clone(&self) -> Self {
+        Self {
+            qr: self.qr.clone(),
+            diag: self.diag.clone(),
+        }
+    }
+}
+
+impl<T: fmt::Debug, R: DimMin<C>, C: Dim> fmt::Debug for QR<T, R, C>
+where
+    DefaultAllocator: Allocator<T, R, C> + Allocator<T, DimMinimum<R, C>>,
+    Owned<T, R, C>: fmt::Debug,
+    Owned<T, DimMinimum<R, C>>: 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<T: ComplexField, R: DimMin<C>, C: Dim> QR<T, R, C>
 where
     DefaultAllocator: Allocator<T, R, C> + Allocator<T, R> + Allocator<T, DimMinimum<R, C>>,
@@ -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<T, D, D>,
          OMatrix<T, D, D>: Deserialize<'de>"))
 )]
-#[derive(Clone, Debug)]
-pub struct Schur<T: ComplexField, D: Dim>
+pub struct Schur<T, D: Dim>
 where
     DefaultAllocator: Allocator<T, D, D>,
 {
@@ -41,13 +42,39 @@ where
     t: OMatrix<T, D, D>,
 }
 
-impl<T: ComplexField, D: Dim> Copy for Schur<T, D>
+impl<T: Copy, D: Dim> Copy for Schur<T, D>
 where
     DefaultAllocator: Allocator<T, D, D>,
-    OMatrix<T, D, D>: Copy,
+    Owned<T, D, D>: Copy,
 {
 }
 
+impl<T: Clone, D: Dim> Clone for Schur<T, D>
+where
+    DefaultAllocator: Allocator<T, D, D>,
+    Owned<T, D, D>: Clone,
+{
+    fn clone(&self) -> Self {
+        Self {
+            q: self.q.clone(),
+            t: self.t.clone(),
+        }
+    }
+}
+
+impl<T: fmt::Debug, D: Dim> fmt::Debug for Schur<T, D>
+where
+    DefaultAllocator: Allocator<T, D, D>,
+    Owned<T, D, D>: 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<T: ComplexField, D: Dim> Schur<T, D>
 where
     D: DimSub<U1>, // 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<T, D, D>, eps: T::RealField, max_niter: usize) -> Option<Self> {
-        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<T, D, D>,
-        work: &mut OVector<T, D>,
+        work: &mut OVector<MaybeUninit<T>, 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<T, D, D>, out: &mut OVector<T, D>) -> bool {
+    fn do_eigenvalues(t: &OMatrix<T, D, D>, out: &mut OVector<MaybeUninit<T>, 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<T, D, D>,
-        out: &mut OVector<NumComplex<T>, D>,
+        out: &mut OVector<MaybeUninit<NumComplex<T>>, D>,
     ) where
         T: RealField,
         DefaultAllocator: Allocator<NumComplex<T>, 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<OVector<T, 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>);
         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<NumComplex<T>, 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<OVector<T, D>> {
+        // 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<NumComplex<T>, 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<T::RealField, DimMinimum<R, C>>: Deserialize<'de>"
     ))
 )]
-#[derive(Clone, Debug)]
 pub struct SVD<T: ComplexField, R: DimMin<C>, C: Dim>
 where
     DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>
@@ -59,12 +60,48 @@ where
     DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>
         + Allocator<T, R, DimMinimum<R, C>>
         + Allocator<T::RealField, DimMinimum<R, C>>,
-    OMatrix<T, R, DimMinimum<R, C>>: Copy,
-    OMatrix<T, DimMinimum<R, C>, C>: Copy,
-    OVector<T::RealField, DimMinimum<R, C>>: Copy,
+    Owned<T, R, DimMinimum<R, C>>: Copy,
+    Owned<T, DimMinimum<R, C>, C>: Copy,
+    Owned<T::RealField, DimMinimum<R, C>>: Copy,
 {
 }
 
+impl<T: ComplexField, R: DimMin<C>, C: Dim> Clone for SVD<T, R, C>
+where
+    DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>
+        + Allocator<T, R, DimMinimum<R, C>>
+        + Allocator<T::RealField, DimMinimum<R, C>>,
+    Owned<T, R, DimMinimum<R, C>>: Clone,
+    Owned<T, DimMinimum<R, C>, C>: Clone,
+    Owned<T::RealField, DimMinimum<R, C>>: Clone,
+{
+    fn clone(&self) -> Self {
+        Self {
+            u: self.u.clone(),
+            v_t: self.v_t.clone(),
+            singular_values: self.singular_values.clone(),
+        }
+    }
+}
+
+impl<T: ComplexField, R: DimMin<C>, C: Dim> fmt::Debug for SVD<T, R, C>
+where
+    DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>
+        + Allocator<T, R, DimMinimum<R, C>>
+        + Allocator<T::RealField, DimMinimum<R, C>>,
+    Owned<T, R, DimMinimum<R, C>>: fmt::Debug,
+    Owned<T, DimMinimum<R, C>, C>: fmt::Debug,
+    Owned<T::RealField, DimMinimum<R, C>>: 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<T: ComplexField, R: DimMin<C>, C: Dim> SVD<T, R, C>
 where
     DimMinimum<R, C>: DimSub<U1>, // 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<T::RealField, D>: Deserialize<'de>,
          OMatrix<T, D, D>: Deserialize<'de>"))
 )]
-#[derive(Clone, Debug)]
 pub struct SymmetricEigen<T: ComplexField, D: Dim>
 where
     DefaultAllocator: Allocator<T, D, D> + Allocator<T::RealField, D>,
@@ -44,11 +45,39 @@ where
 impl<T: ComplexField, D: Dim> Copy for SymmetricEigen<T, D>
 where
     DefaultAllocator: Allocator<T, D, D> + Allocator<T::RealField, D>,
-    OMatrix<T, D, D>: Copy,
-    OVector<T::RealField, D>: Copy,
+    Owned<T, D, D>: Copy,
+    Owned<T::RealField, D>: Copy,
 {
 }
 
+impl<T: ComplexField, D: Dim> Clone for SymmetricEigen<T, D>
+where
+    DefaultAllocator: Allocator<T, D, D> + Allocator<T::RealField, D>,
+    Owned<T, D, D>: Clone,
+    Owned<T::RealField, D>: Clone,
+{
+    fn clone(&self) -> Self {
+        Self {
+            eigenvectors: self.eigenvectors.clone(),
+            eigenvalues: self.eigenvalues.clone(),
+        }
+    }
+}
+
+impl<T: ComplexField, D: Dim> fmt::Debug for SymmetricEigen<T, D>
+where
+    DefaultAllocator: Allocator<T, D, D> + Allocator<T::RealField, D>,
+    Owned<T, D, D>: fmt::Debug,
+    Owned<T::RealField, D>: 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<T: ComplexField, D: Dim> SymmetricEigen<T, D>
 where
     DefaultAllocator: Allocator<T, D, D> + Allocator<T::RealField, D>,
@@ -270,7 +299,10 @@ where
     ///
     /// This is useful if some of the eigenvalues have been manually modified.
     #[must_use]
-    pub fn recompose(&self) -> OMatrix<T, D, D> {
+    pub fn recompose(&self) -> OMatrix<T, D, D>
+    where
+        Owned<T, D, D>: 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<T, D, D>: Deserialize<'de>,
          OVector<T, DimDiff<D, U1>>: Deserialize<'de>"))
 )]
-#[derive(Clone, Debug)]
-pub struct SymmetricTridiagonal<T: ComplexField, D: DimSub<U1>>
+pub struct SymmetricTridiagonal<T, D: DimSub<U1>>
 where
     DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
 {
@@ -34,14 +36,42 @@ where
     off_diagonal: OVector<T, DimDiff<D, U1>>,
 }
 
-impl<T: ComplexField, D: DimSub<U1>> Copy for SymmetricTridiagonal<T, D>
+impl<T: Copy, D: DimSub<U1>> Copy for SymmetricTridiagonal<T, D>
 where
     DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
-    OMatrix<T, D, D>: Copy,
-    OVector<T, DimDiff<D, U1>>: Copy,
+    Owned<T, D, D>: Copy,
+    Owned<T, DimDiff<D, U1>>: Copy,
 {
 }
 
+impl<T: Clone, D: DimSub<U1>> Clone for SymmetricTridiagonal<T, D>
+where
+    DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
+    Owned<T, D, D>: Clone,
+    Owned<T, DimDiff<D, U1>>: Clone,
+{
+    fn clone(&self) -> Self {
+        Self {
+            tri: self.tri.clone(),
+            off_diagonal: self.off_diagonal.clone(),
+        }
+    }
+}
+
+impl<T: fmt::Debug, D: DimSub<U1>> fmt::Debug for SymmetricTridiagonal<T, D>
+where
+    DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
+    Owned<T, D, D>: fmt::Debug,
+    Owned<T, DimDiff<D, U1>>: 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<T: ComplexField, D: DimSub<U1>> SymmetricTridiagonal<T, D>
 where
     DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
@@ -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<T, D>: Deserialize<'de>, OMatrix<T, D, D>: Deserialize<'de>"
     ))
 )]
-#[derive(Clone, Debug)]
-pub struct UDU<T: RealField, D: Dim>
+pub struct UDU<T, D: Dim>
 where
     DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
 {
@@ -30,14 +31,42 @@ where
     pub d: OVector<T, D>,
 }
 
-impl<T: RealField, D: Dim> Copy for UDU<T, D>
+impl<T: Copy, D: Dim> Copy for UDU<T, D>
 where
     DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
-    OVector<T, D>: Copy,
-    OMatrix<T, D, D>: Copy,
+    Owned<T, D>: Copy,
+    Owned<T, D, D>: Copy,
 {
 }
 
+impl<T: Clone, D: Dim> Clone for UDU<T, D>
+where
+    DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
+    Owned<T, D>: Clone,
+    Owned<T, D, D>: Clone,
+{
+    fn clone(&self) -> Self {
+        Self {
+            u: self.u.clone(),
+            d: self.d.clone(),
+        }
+    }
+}
+
+impl<T: fmt::Debug, D: Dim> fmt::Debug for UDU<T, D>
+where
+    DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
+    Owned<T, D>: fmt::Debug,
+    Owned<T, D, D>: 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<T: RealField, D: Dim> UDU<T, D>
 where
     DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,