diff --git a/src/base/default_allocator.rs b/src/base/default_allocator.rs
index 0cd6874b..519f85f3 100644
--- a/src/base/default_allocator.rs
+++ b/src/base/default_allocator.rs
@@ -31,6 +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
 /// An allocator based on `GenericArray` and `VecStorage` for statically-sized and dynamically-sized
 /// matrices respectively.
 pub struct DefaultAllocator;
diff --git a/src/base/matrix.rs b/src/base/matrix.rs
index 71c3b38e..d13a467e 100644
--- a/src/base/matrix.rs
+++ b/src/base/matrix.rs
@@ -152,7 +152,7 @@ pub type MatrixCross<T, R1, C1, R2, C2> =
 /// dynamically-sized column vector should be represented as a `Matrix<T, Dynamic, U1, S>` (given
 /// some concrete types for `T` and a compatible data storage type `S`).
 #[repr(C)]
-#[derive(Clone, Copy)]
+#[derive(Clone, Copy, Debug)]
 pub struct Matrix<T, R, C, S> {
     /// The data storage that contains all the matrix components. Disappointed?
     ///
@@ -192,15 +192,6 @@ pub struct Matrix<T, R, C, S> {
     _phantoms: PhantomData<(T, R, C)>,
 }
 
-impl<T, R: Dim, C: Dim, S: fmt::Debug> fmt::Debug for Matrix<T, R, C, S> {
-    fn fmt(&self, formatter: &mut fmt::Formatter) -> Result<(), fmt::Error> {
-        formatter
-            .debug_struct("Matrix")
-            .field("data", &self.data)
-            .finish()
-    }
-}
-
 impl<T, R: Dim, C: Dim, S: Default> Default for Matrix<T, R, C, S> {
     fn default() -> Self {
         Matrix {
diff --git a/src/base/unit.rs b/src/base/unit.rs
index 8346d2ed..f656b247 100644
--- a/src/base/unit.rs
+++ b/src/base/unit.rs
@@ -228,7 +228,7 @@ impl<T> Unit<T> {
     /// Wraps the given reference, assuming it is already normalized.
     #[inline]
     pub fn from_ref_unchecked(value: &T) -> &Self {
-        unsafe { &*(value as *const T as *const Self) }
+        unsafe { &*(value as *const _ as *const Self) }
     }
 
     /// Retrieves the underlying value.
@@ -331,7 +331,7 @@ impl<T> Deref for Unit<T> {
 
     #[inline]
     fn deref(&self) -> &T {
-        unsafe { &*(self as *const Self as *const T) }
+        unsafe { &*(self as *const _ as *const T) }
     }
 }
 
diff --git a/src/geometry/dual_quaternion_ops.rs b/src/geometry/dual_quaternion_ops.rs
index 7d07ec2c..4f1e58e3 100644
--- a/src/geometry/dual_quaternion_ops.rs
+++ b/src/geometry/dual_quaternion_ops.rs
@@ -59,14 +59,14 @@ use std::ops::{
 impl<T> AsRef<[T; 8]> for DualQuaternion<T> {
     #[inline]
     fn as_ref(&self) -> &[T; 8] {
-        unsafe { &*(self as *const Self as *const [T; 8]) }
+        unsafe { &*(self as *const _ as *const [T; 8]) }
     }
 }
 
 impl<T> AsMut<[T; 8]> for DualQuaternion<T> {
     #[inline]
     fn as_mut(&mut self) -> &mut [T; 8] {
-        unsafe { &mut *(self as *mut Self as *mut [T; 8]) }
+        unsafe { &mut *(self as *mut _ as *mut [T; 8]) }
     }
 }
 
diff --git a/src/geometry/orthographic.rs b/src/geometry/orthographic.rs
index e9546cdd..98fd6b0d 100644
--- a/src/geometry/orthographic.rs
+++ b/src/geometry/orthographic.rs
@@ -18,26 +18,27 @@ use crate::base::{Matrix4, Vector, Vector3};
 use crate::geometry::{Point3, Projective3};
 
 /// A 3D orthographic projection stored as a homogeneous 4x4 matrix.
+#[repr(C)]
 pub struct Orthographic3<T> {
     matrix: Matrix4<T>,
 }
 
-impl<T: RealField> Copy for Orthographic3<T> {}
+impl<T: Copy> Copy for Orthographic3<T> {}
 
-impl<T: RealField> Clone for Orthographic3<T> {
+impl<T: Clone> Clone for Orthographic3<T> {
     #[inline]
     fn clone(&self) -> Self {
         Self::from_matrix_unchecked(self.matrix)
     }
 }
 
-impl<T: RealField> fmt::Debug for Orthographic3<T> {
-    fn fmt(&self, f: &mut fmt::Formatter) -> Result<(), fmt::Error> {
+impl<T: fmt::Debug> fmt::Debug for Orthographic3<T> {
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
         self.matrix.fmt(f)
     }
 }
 
-impl<T: RealField> PartialEq for Orthographic3<T> {
+impl<T: PartialEq> PartialEq for Orthographic3<T> {
     #[inline]
     fn eq(&self, right: &Self) -> bool {
         self.matrix == right.matrix
@@ -45,7 +46,7 @@ impl<T: RealField> PartialEq for Orthographic3<T> {
 }
 
 #[cfg(feature = "serde-serialize-no-std")]
-impl<T: RealField + Serialize> Serialize for Orthographic3<T> {
+impl<T: Serialize> Serialize for Orthographic3<T> {
     fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
     where
         S: Serializer,
@@ -55,7 +56,7 @@ impl<T: RealField + Serialize> Serialize for Orthographic3<T> {
 }
 
 #[cfg(feature = "serde-serialize-no-std")]
-impl<'a, T: RealField + Deserialize<'a>> Deserialize<'a> for Orthographic3<T> {
+impl<'a, T: Deserialize<'a>> Deserialize<'a> for Orthographic3<T> {
     fn deserialize<Des>(deserializer: Des) -> Result<Self, Des::Error>
     where
         Des: Deserializer<'a>,
@@ -66,7 +67,8 @@ impl<'a, T: RealField + Deserialize<'a>> Deserialize<'a> for Orthographic3<T> {
     }
 }
 
-impl<T: RealField> Orthographic3<T> {
+/// # Basic methods and casts.
+impl<T> Orthographic3<T> {
     /// Creates a new orthographic projection matrix.
     ///
     /// This follows the OpenGL convention, so this will flip the `z` axis.
@@ -110,8 +112,11 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_relative_eq!(proj.project_point(&p8), Point3::new(-1.0, -1.0, -1.0));
     /// ```
     #[inline]
-    pub fn new(left: T, right: T, bottom: T, top: T, znear: T, zfar: T) -> Self {
-        let matrix = Matrix4::<T>::identity();
+    pub fn new(left: T, right: T, bottom: T, top: T, znear: T, zfar: T) -> Self
+    where
+        T: RealField,
+    {
+        let matrix = Matrix4::identity();
         let mut res = Self::from_matrix_unchecked(matrix);
 
         res.set_left_and_right(left, right);
@@ -145,7 +150,10 @@ impl<T: RealField> Orthographic3<T> {
 
     /// Creates a new orthographic projection matrix from an aspect ratio and the vertical field of view.
     #[inline]
-    pub fn from_fov(aspect: T, vfov: T, znear: T, zfar: T) -> Self {
+    pub fn from_fov(aspect: T, vfov: T, znear: T, zfar: T) -> Self
+    where
+        T: RealField,
+    {
         assert!(
             znear != zfar,
             "The far plane must not be equal to the near plane."
@@ -188,7 +196,10 @@ impl<T: RealField> Orthographic3<T> {
     /// ```
     #[inline]
     #[must_use]
-    pub fn inverse(&self) -> Matrix4<T> {
+    pub fn inverse(&self) -> Matrix4<T>
+    where
+        T: RealField,
+    {
         let mut res = self.to_homogeneous();
 
         let inv_m11 = T::one() / self.matrix[(0, 0)];
@@ -257,7 +268,8 @@ impl<T: RealField> Orthographic3<T> {
     #[inline]
     #[must_use]
     pub fn as_projective(&self) -> &Projective3<T> {
-        unsafe { &*(self as *const Orthographic3<T> as *const Projective3<T>) }
+        // Safety: Self and Projective3 are both #[repr(C)] of a matrix.
+        unsafe { &*(self as *const _ as *const Projective3<T>) }
     }
 
     /// This transformation seen as a `Projective3`.
@@ -301,7 +313,10 @@ impl<T: RealField> Orthographic3<T> {
     pub fn unwrap(self) -> Matrix4<T> {
         self.matrix
     }
+}
 
+/// # Mathematical methods.
+impl<T: RealField> Orthographic3<T> {
     /// The left offset of the view cuboid.
     ///
     /// ```
diff --git a/src/geometry/perspective.rs b/src/geometry/perspective.rs
index ba8368a2..73023080 100644
--- a/src/geometry/perspective.rs
+++ b/src/geometry/perspective.rs
@@ -33,7 +33,7 @@ impl<T: RealField> Clone for Perspective3<T> {
 }
 
 impl<T: RealField> fmt::Debug for Perspective3<T> {
-    fn fmt(&self, f: &mut fmt::Formatter) -> Result<(), fmt::Error> {
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
         self.matrix.fmt(f)
     }
 }
@@ -139,7 +139,7 @@ impl<T: RealField> Perspective3<T> {
     #[inline]
     #[must_use]
     pub fn as_projective(&self) -> &Projective3<T> {
-        unsafe { &*(self as *const Perspective3<T> as *const Projective3<T>) }
+        unsafe { &*(self as *const _ as *const Projective3<T>) }
     }
 
     /// This transformation seen as a `Projective3`.
diff --git a/src/geometry/quaternion_coordinates.rs b/src/geometry/quaternion_coordinates.rs
index cb16e59e..ba887f63 100644
--- a/src/geometry/quaternion_coordinates.rs
+++ b/src/geometry/quaternion_coordinates.rs
@@ -12,13 +12,13 @@ impl<T: Scalar + SimdValue> Deref for Quaternion<T> {
 
     #[inline]
     fn deref(&self) -> &Self::Target {
-        unsafe { &*(self as *const Self as *const Self::Target) }
+        unsafe { &*(self as *const _ as *const Self::Target) }
     }
 }
 
 impl<T: Scalar + SimdValue> DerefMut for Quaternion<T> {
     #[inline]
     fn deref_mut(&mut self) -> &mut Self::Target {
-        unsafe { &mut *(self as *mut Self as *mut Self::Target) }
+        unsafe { &mut *(self as *mut _ as *mut Self::Target) }
     }
 }
diff --git a/src/geometry/translation_coordinates.rs b/src/geometry/translation_coordinates.rs
index 80267e06..44a4c8f2 100644
--- a/src/geometry/translation_coordinates.rs
+++ b/src/geometry/translation_coordinates.rs
@@ -18,14 +18,14 @@ macro_rules! deref_impl(
 
             #[inline]
             fn deref(&self) -> &Self::Target {
-                unsafe { &*(self as *const Translation<T, $D> as *const Self::Target) }
+                unsafe { &*(self as *const _ as *const Self::Target) }
             }
         }
 
         impl<T: Scalar> DerefMut for Translation<T, $D> {
             #[inline]
             fn deref_mut(&mut self) -> &mut Self::Target {
-                unsafe { &mut *(self as *mut Translation<T, $D> as *mut Self::Target) }
+                unsafe { &mut *(self as *mut _ as *mut Self::Target) }
             }
         }
     }
diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs
index ac40331f..46bb9029 100644
--- a/src/linalg/bidiagonal.rs
+++ b/src/linalg/bidiagonal.rs
@@ -7,7 +7,7 @@ use crate::allocator::Allocator;
 use crate::base::{DefaultAllocator, Matrix, OMatrix, OVector, Unit};
 use crate::dimension::{Const, Dim, DimDiff, DimMin, DimMinimum, DimSub, U1};
 use crate::storage::{Owned, Storage};
-use crate::{Dynamic, };
+use crate::Dynamic;
 use simba::scalar::ComplexField;
 
 use crate::geometry::Reflection;
diff --git a/src/linalg/col_piv_qr.rs b/src/linalg/col_piv_qr.rs
index 1a56d2cb..1d01f294 100644
--- a/src/linalg/col_piv_qr.rs
+++ b/src/linalg/col_piv_qr.rs
@@ -30,7 +30,6 @@ use crate::linalg::{householder, PermutationSequence};
          PermutationSequence<DimMinimum<R, C>>: Deserialize<'de>,
          OVector<T, DimMinimum<R, C>>: Deserialize<'de>"))
 )]
-#[derive(Clone, Debug)]
 pub struct ColPivQR<T: ComplexField, R: DimMin<C>, C: Dim>
 where
     DefaultAllocator: Allocator<T, R, C>
@@ -53,6 +52,24 @@ where
 {
 }
 
+impl<T: ComplexField, R: DimMin<C>, C: Dim> Clone for ColPivQR<T, R, C>
+where
+    DefaultAllocator: Allocator<T, R, C>
+        + Allocator<T, DimMinimum<R, C>>
+        + Allocator<(usize, usize), DimMinimum<R, C>>,
+    OMatrix<T, R, C>: Clone,
+    PermutationSequence<DimMinimum<R, C>>: Clone,
+    OVector<T, DimMinimum<R, C>>: Clone,
+{
+    fn clone(&self) -> Self {
+        Self {
+            col_piv_qr: self.col_piv_qr.clone(),
+            p: self.p.clone(),
+            diag: self.diag.clone(),
+        }
+    }
+}
+
 impl<T: ComplexField, R: DimMin<C>, C: Dim> ColPivQR<T, R, C>
 where
     DefaultAllocator: Allocator<T, R, C>
@@ -66,14 +83,13 @@ where
         let min_nrows_ncols = nrows.min(ncols);
         let mut p = PermutationSequence::identity_generic(min_nrows_ncols);
 
-        let mut diag =
-            unsafe { crate::unimplemented_or_uninitialized_generic!(min_nrows_ncols, Const::<1>) };
+        let mut diag = Matrix::new_uninitialized_generic(min_nrows_ncols, Const::<1>);
 
         if min_nrows_ncols.value() == 0 {
             return ColPivQR {
                 col_piv_qr: matrix,
                 p,
-                diag,
+                diag: unsafe { diag.assume_init() },
             };
         }
 
@@ -83,13 +99,13 @@ where
             matrix.swap_columns(i, col_piv);
             p.append_permutation(i, col_piv);
 
-            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);
         }
 
         ColPivQR {
             col_piv_qr: matrix,
             p,
-            diag,
+            diag:unsafe{diag.assume_init()},
         }
     }
 
diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index c2816ff0..4fc5b460 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -1,14 +1,11 @@
 //! This module provides the matrix exponent (exp) function to square matrices.
 //!
-use crate::{
-    base::{
+use crate::{ComplexField, OMatrix, RealField, base::{
         allocator::Allocator,
         dimension::{Const, Dim, DimMin, DimMinimum},
         storage::Storage,
         DefaultAllocator,
-    },
-    convert, try_convert, ComplexField, OMatrix, RealField,
-};
+    }, convert, storage::Owned, try_convert};
 
 use crate::num::Zero;
 
@@ -433,6 +430,7 @@ where
         + Allocator<T, D>
         + Allocator<T::RealField, D>
         + Allocator<T::RealField, D, D>,
+    Owned<T, D, D>: Clone,
 {
     /// Computes exponential of this matrix
     #[must_use]
diff --git a/src/linalg/full_piv_lu.rs b/src/linalg/full_piv_lu.rs
index f08af55c..71e0755e 100644
--- a/src/linalg/full_piv_lu.rs
+++ b/src/linalg/full_piv_lu.rs
@@ -1,3 +1,5 @@
+use std::fmt;
+
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Serialize};
 
@@ -27,8 +29,7 @@ use crate::linalg::PermutationSequence;
          OMatrix<T, R, C>: Deserialize<'de>,
          PermutationSequence<DimMinimum<R, C>>: Deserialize<'de>"))
 )]
-#[derive(Clone, Debug)]
-pub struct FullPivLU<T: ComplexField, R: DimMin<C>, C: Dim>
+pub struct FullPivLU<T, R: DimMin<C>, C: Dim>
 where
     DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
 {
@@ -40,11 +41,41 @@ where
 impl<T: ComplexField, R: DimMin<C>, C: Dim> Copy for FullPivLU<T, R, C>
 where
     DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
-    OMatrix<T, R, C>: Copy,
     PermutationSequence<DimMinimum<R, C>>: Copy,
+    OMatrix<T, R, C>: Copy,
 {
 }
 
+impl<T: ComplexField, R: DimMin<C>, C: Dim> Clone for FullPivLU<T, R, C>
+where
+    DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
+    PermutationSequence<DimMinimum<R, C>>: Clone,
+    OMatrix<T, R, C>: Clone,
+{
+    fn clone(&self) -> Self {
+        Self {
+            lu: self.lu.clone(),
+            p: self.p.clone(),
+            q: self.q.clone(),
+        }
+    }
+}
+
+impl<T: ComplexField, R: DimMin<C>, C: Dim> fmt::Debug for FullPivLU<T, R, C>
+where
+    DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
+    PermutationSequence<DimMinimum<R, C>>: fmt::Debug,
+    OMatrix<T, R, C>: fmt::Debug,
+{
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        f.debug_struct("FullPivLU")
+            .field("lu", &self.lu)
+            .field("p", &self.p)
+            .field("q", &self.q)
+            .finish()
+    }
+}
+
 impl<T: ComplexField, R: DimMin<C>, C: Dim> FullPivLU<T, R, C>
 where
     DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
diff --git a/src/linalg/hessenberg.rs b/src/linalg/hessenberg.rs
index 6b8ecfee..6a4260bf 100644
--- a/src/linalg/hessenberg.rs
+++ b/src/linalg/hessenberg.rs
@@ -1,10 +1,14 @@
+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 crate::Matrix;
 use simba::scalar::ComplexField;
 
 use crate::linalg::householder;
@@ -25,7 +29,6 @@ use crate::linalg::householder;
          OMatrix<T, D, D>: Deserialize<'de>,
          OVector<T, DimDiff<D, U1>>: Deserialize<'de>"))
 )]
-#[derive(Clone, Debug)]
 pub struct Hessenberg<T: ComplexField, D: DimSub<U1>>
 where
     DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
@@ -37,20 +40,46 @@ where
 impl<T: ComplexField, D: DimSub<U1>> Copy for Hessenberg<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: ComplexField, D: DimSub<U1>> Clone for Hessenberg<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 {
+            hess: self.hess.clone(),
+            subdiag: self.subdiag.clone(),
+        }
+    }
+}
+
+impl<T: ComplexField, D: DimSub<U1>> fmt::Debug for Hessenberg<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("Hessenberg")
+            .field("hess", &self.hess)
+            .field("subdiag", &self.subdiag)
+            .finish()
+    }
+}
+
 impl<T: ComplexField, D: DimSub<U1>> Hessenberg<T, D>
 where
     DefaultAllocator: Allocator<T, D, D> + Allocator<T, D> + Allocator<T, DimDiff<D, U1>>,
 {
     /// Computes the Hessenberg decomposition using householder reflections.
     pub fn new(hess: OMatrix<T, D, D>) -> Self {
-        let mut work = unsafe {
-            crate::unimplemented_or_uninitialized_generic!(hess.data.shape().0, Const::<1>)
-        };
+        let mut work = OVector::new_uninitialized_generic(hess.data.shape().0, Const::<1>);
         Self::new_with_workspace(hess, &mut work)
     }
 
@@ -58,7 +87,10 @@ where
     ///
     /// The workspace containing `D` elements must be provided but its content does not have to be
     /// initialized.
-    pub fn new_with_workspace(mut hess: OMatrix<T, D, D>, work: &mut OVector<T, D>) -> Self {
+    pub fn new_with_workspace(
+        mut hess: OMatrix<T, D, D>,
+        work: &mut OVector<MaybeUninit<T>, D>,
+    ) -> Self {
         assert!(
             hess.is_square(),
             "Cannot compute the hessenberg decomposition of a non-square matrix."
@@ -76,19 +108,29 @@ where
             "Hessenberg: invalid workspace size."
         );
 
-        let mut subdiag = unsafe {
-            crate::unimplemented_or_uninitialized_generic!(dim.sub(Const::<1>), Const::<1>)
-        };
+        let mut subdiag = Matrix::new_uninitialized_generic(dim.sub(Const::<1>), Const::<1>);
 
         if dim.value() == 0 {
-            return Hessenberg { hess, subdiag };
+            return Self {
+                hess,
+                subdiag: unsafe { subdiag.assume_init() },
+            };
         }
 
         for ite in 0..dim.value() - 1 {
-            householder::clear_column_unchecked(&mut hess, &mut subdiag[ite], ite, 1, Some(work));
+            householder::clear_column_unchecked(
+                &mut hess,
+                subdiag[ite].as_mut_ptr(),
+                ite,
+                1,
+                Some(work),
+            );
         }
 
-        Hessenberg { hess, subdiag }
+        Self {
+            hess,
+            subdiag: unsafe { subdiag.assume_init() },
+        }
     }
 
     /// Retrieves `(q, h)` with `q` the orthogonal matrix of this decomposition and `h` the
@@ -117,7 +159,10 @@ where
     /// This is less efficient than `.unpack_h()` as it allocates a new matrix.
     #[inline]
     #[must_use]
-    pub fn h(&self) -> OMatrix<T, D, D> {
+    pub fn h(&self) -> OMatrix<T, D, D>
+    where
+        Owned<T, D, D>: Clone,
+    {
         let dim = self.hess.nrows();
         let mut res = self.hess.clone();
         res.fill_lower_triangle(T::zero(), 2);
diff --git a/src/linalg/householder.rs b/src/linalg/householder.rs
index c53bc4b4..cb65900a 100644
--- a/src/linalg/householder.rs
+++ b/src/linalg/householder.rs
@@ -51,7 +51,7 @@ pub fn clear_column_unchecked<T: ComplexField, R: Dim, C: Dim>(
     diag_elt: *mut T,
     icol: usize,
     shift: usize,
-    bilateral: Option<&mut OVector<T, R>>,
+    bilateral: Option<&mut OVector<MaybeUninit<T>, R>>,
 ) where
     DefaultAllocator: Allocator<T, R, C> + Allocator<T, R>,
 {
@@ -88,11 +88,14 @@ pub fn clear_row_unchecked<T: ComplexField, R: Dim, C: Dim>(
 {
     let (mut top, mut bottom) = matrix.rows_range_pair_mut(irow, irow + 1..);
     let mut axis = axis_packed.rows_range_mut(irow + shift..);
-    axis.tr_copy_from(&top.columns_range(irow + shift..));
+    axis.tr_copy_init_from(&top.columns_range(irow + shift..));
+    let mut axis = unsafe { axis.assume_init_mut() };
 
     let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis);
     axis.conjugate_mut(); // So that reflect_rows actually cancels the first row.
-  unsafe{  *diag_elt = reflection_norm;}
+    unsafe {
+        *diag_elt = reflection_norm;
+    }
 
     if not_zero {
         let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
diff --git a/src/linalg/lu.rs b/src/linalg/lu.rs
index 36a00807..8b4fb7c3 100644
--- a/src/linalg/lu.rs
+++ b/src/linalg/lu.rs
@@ -1,3 +1,6 @@
+use std::fmt;
+use std::mem;
+
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Serialize};
 
@@ -5,9 +8,8 @@ use crate::allocator::{Allocator, Reallocator};
 use crate::base::{DefaultAllocator, Matrix, OMatrix, Scalar};
 use crate::constraint::{SameNumberOfRows, ShapeConstraint};
 use crate::dimension::{Dim, DimMin, DimMinimum};
-use crate::storage::{Storage, StorageMut};
+use crate::storage::{Owned, Storage, StorageMut};
 use simba::scalar::{ComplexField, Field};
-use std::mem;
 
 use crate::linalg::PermutationSequence;
 
@@ -27,8 +29,7 @@ use crate::linalg::PermutationSequence;
          OMatrix<T, R, C>: Deserialize<'de>,
          PermutationSequence<DimMinimum<R, C>>: Deserialize<'de>"))
 )]
-#[derive(Clone, Debug)]
-pub struct LU<T: ComplexField, R: DimMin<C>, C: Dim>
+pub struct LU<T, R: DimMin<C>, C: Dim>
 where
     DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
 {
@@ -36,14 +37,42 @@ where
     p: PermutationSequence<DimMinimum<R, C>>,
 }
 
-impl<T: ComplexField, R: DimMin<C>, C: Dim> Copy for LU<T, R, C>
+impl<T: Copy, R: DimMin<C>, C: Dim> Copy for LU<T, R, C>
 where
     DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
-    OMatrix<T, R, C>: Copy,
     PermutationSequence<DimMinimum<R, C>>: Copy,
+    Owned<T, R, C>: Copy,
 {
 }
 
+impl<T: Clone, R: DimMin<C>, C: Dim> Clone for LU<T, R, C>
+where
+    DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
+    PermutationSequence<DimMinimum<R, C>>: Clone,
+    Owned<T, R, C>: Clone,
+{
+    fn clone(&self) -> Self {
+        Self {
+            lu: self.lu.clone(),
+            p: self.p.clone(),
+        }
+    }
+}
+
+impl<T: ComplexField, R: DimMin<C>, C: Dim> fmt::Debug for LU<T, R, C>
+where
+    DefaultAllocator: Allocator<T, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
+    PermutationSequence<DimMinimum<R, C>>: fmt::Debug,
+    Owned<T, R, C>: fmt::Debug,
+{
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        f.debug_struct("LU")
+            .field("lu", &self.lu)
+            .field("p", &self.p)
+            .finish()
+    }
+}
+
 /// Performs a LU decomposition to overwrite `out` with the inverse of `matrix`.
 ///
 /// If `matrix` is not invertible, `false` is returned and `out` may contain invalid data.
diff --git a/src/linalg/permutation_sequence.rs b/src/linalg/permutation_sequence.rs
index a088c458..e4594520 100644
--- a/src/linalg/permutation_sequence.rs
+++ b/src/linalg/permutation_sequence.rs
@@ -1,3 +1,4 @@
+use std::fmt;
 use std::mem::MaybeUninit;
 
 #[cfg(feature = "serde-serialize-no-std")]
@@ -10,8 +11,10 @@ use crate::allocator::Allocator;
 use crate::base::{DefaultAllocator, Matrix, OVector, Scalar};
 #[cfg(any(feature = "std", feature = "alloc"))]
 use crate::dimension::Dynamic;
-use crate::dimension::{ Dim, DimName};
-use crate::storage::StorageMut;
+use crate::dimension::{Dim, DimName};
+use crate::iter::MatrixIter;
+use crate::storage::{Owned, StorageMut};
+use crate::{Const, U1};
 
 /// A sequence of row or column permutations.
 #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))]
@@ -25,7 +28,6 @@ use crate::storage::StorageMut;
     serde(bound(deserialize = "DefaultAllocator: Allocator<(usize, usize), D>,
          OVector<(usize, usize), D>: Deserialize<'de>"))
 )]
-#[derive(Clone, Debug)]
 pub struct PermutationSequence<D: Dim>
 where
     DefaultAllocator: Allocator<(usize, usize), D>,
@@ -41,6 +43,32 @@ where
 {
 }
 
+impl<D: Dim> Clone for PermutationSequence<D>
+where
+    DefaultAllocator: Allocator<(usize, usize), D>,
+    OVector<MaybeUninit<(usize, usize)>, D>: Clone,
+{
+    fn clone(&self) -> Self {
+        Self {
+            len: self.len,
+            ipiv: self.ipiv.clone(),
+        }
+    }
+}
+
+impl<D: Dim> fmt::Debug for PermutationSequence<D>
+where
+    DefaultAllocator: Allocator<(usize, usize), D>,
+    OVector<MaybeUninit<(usize, usize)>, D>: fmt::Debug,
+{
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        f.debug_struct("PermutationSequence")
+            .field("len", &self.len)
+            .field("ipiv", &self.ipiv)
+            .finish()
+    }
+}
+
 impl<D: DimName> PermutationSequence<D>
 where
     DefaultAllocator: Allocator<(usize, usize), D>,
@@ -74,7 +102,7 @@ where
         unsafe {
             Self {
                 len: 0,
-                ipiv: OVector::new_uninitialized(dim),
+                ipiv: OVector::new_uninitialized_generic(dim, Const::<1>),
             }
         }
     }
@@ -88,7 +116,7 @@ where
                 self.len < self.ipiv.len(),
                 "Maximum number of permutations exceeded."
             );
-            self.ipiv[self.len] = (i, i2);
+            self.ipiv[self.len] = MaybeUninit::new((i, i2));
             self.len += 1;
         }
     }
@@ -99,8 +127,8 @@ where
     where
         S2: StorageMut<T, R2, C2>,
     {
-        for i in self.ipiv.rows_range(..self.len).iter().map(MaybeUninit::assume_init) {
-            rhs.swap_rows(i.0, i.1)
+        for perm in self.iter() {
+            rhs.swap_rows(perm.0, perm.1)
         }
     }
 
@@ -110,8 +138,8 @@ where
     where
         S2: StorageMut<T, R2, C2>,
     {
-        for i in 0..self.len {
-            let (i1, i2) = self.ipiv[self.len - i - 1];
+        for perm in self.iter().rev() {
+            let (i1, i2) = perm;
             rhs.swap_rows(i1, i2)
         }
     }
@@ -122,8 +150,8 @@ where
     where
         S2: StorageMut<T, R2, C2>,
     {
-        for i in self.ipiv.rows_range(..self.len).iter() {
-            rhs.swap_columns(i.0, i.1)
+        for perm in self.iter() {
+            rhs.swap_columns(perm.0, perm.1)
         }
     }
 
@@ -135,8 +163,8 @@ where
     ) where
         S2: StorageMut<T, R2, C2>,
     {
-        for i in 0..self.len {
-            let (i1, i2) = self.ipiv[self.len - i - 1];
+        for perm in self.iter().rev() {
+            let (i1, i2) = perm;
             rhs.swap_columns(i1, i2)
         }
     }
@@ -163,4 +191,27 @@ where
             -T::one()
         }
     }
+
+    /// Iterates over the permutations that have been initialized.
+    pub fn iter(
+        &self,
+    ) -> std::iter::Map<
+        std::iter::Copied<
+            std::iter::Take<
+                MatrixIter<
+                    MaybeUninit<(usize, usize)>,
+                    D,
+                    U1,
+                    Owned<MaybeUninit<(usize, usize)>, D, U1>,
+                >,
+            >,
+        >,
+        impl FnMut(MaybeUninit<(usize, usize)>) -> (usize, usize),
+    > {
+        self.ipiv
+            .iter()
+            .take(self.len)
+            .copied()
+            .map(|e| unsafe { e.assume_init() })
+    }
 }