From 10b5dc9bb6e1fd458a5e94c07d665a0a01bb58a5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Violeta=20Hern=C3=A1ndez?= <vi.hdz.p@gmail.com>
Date: Sat, 17 Jul 2021 20:19:20 -0500
Subject: [PATCH] Many miscellaneous improvements throughout

---
 nalgebra-lapack/src/eigen.rs            |  37 ++++++++-
 src/base/allocator.rs                   |   7 +-
 src/base/array_storage.rs               |   8 +-
 src/base/blas.rs                        |  76 +++++++++++------
 src/base/conversion.rs                  |  32 ++++----
 src/base/default_allocator.rs           |   3 +-
 src/base/dimension.rs                   |   9 +-
 src/base/indexing.rs                    |   4 +-
 src/base/matrix.rs                      | 104 ++++++++++++++----------
 src/base/matrix_slice.rs                |  33 ++++----
 src/base/ops.rs                         |  17 ++--
 src/base/scalar.rs                      |  10 ++-
 src/base/unit.rs                        |   2 +-
 src/base/vec_storage.rs                 |   1 -
 src/geometry/dual_quaternion.rs         |   1 +
 src/geometry/dual_quaternion_ops.rs     |   4 +-
 src/geometry/isometry.rs                |   1 -
 src/geometry/orthographic.rs            |   8 +-
 src/geometry/perspective.rs             |   3 +-
 src/geometry/point.rs                   |   2 +-
 src/geometry/point_construction.rs      |   2 +-
 src/geometry/quaternion.rs              |   2 +-
 src/geometry/quaternion_coordinates.rs  |   5 +-
 src/geometry/reflection.rs              |   8 +-
 src/geometry/rotation.rs                |   4 +-
 src/geometry/similarity.rs              |   1 -
 src/geometry/transform.rs               |   2 +-
 src/geometry/translation.rs             |   2 +-
 src/geometry/translation_coordinates.rs |   4 +-
 src/linalg/bidiagonal.rs                |  85 ++++++++++---------
 src/linalg/col_piv_qr.rs                |  27 ++++--
 src/linalg/hessenberg.rs                |  37 +++++----
 src/linalg/householder.rs               |  32 ++++++--
 src/linalg/pow.rs                       |  17 ++--
 src/linalg/qr.rs                        |  14 +++-
 src/proptest/mod.rs                     |   4 +-
 36 files changed, 374 insertions(+), 234 deletions(-)

diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs
index 4347cb03..49fb72b4 100644
--- a/nalgebra-lapack/src/eigen.rs
+++ b/nalgebra-lapack/src/eigen.rs
@@ -1,3 +1,5 @@
+use std::fmt;
+
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Serialize};
 
@@ -32,8 +34,7 @@ use lapack;
          OMatrix<T, D, D>: Deserialize<'de>")
     )
 )]
-#[derive(Clone, Debug)]
-pub struct Eigen<T: Scalar, D: Dim>
+pub struct Eigen<T, D: Dim>
 where
     DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
 {
@@ -45,7 +46,7 @@ where
     pub left_eigenvectors: Option<OMatrix<T, D, D>>,
 }
 
-impl<T: Scalar + Copy, D: Dim> Copy for Eigen<T, D>
+impl<T: Copy, D: Dim> Copy for Eigen<T, D>
 where
     DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
     OVector<T, D>: Copy,
@@ -53,6 +54,36 @@ where
 {
 }
 
+impl<T: Clone, D: Dim> Clone for Eigen<T, D>
+where
+    DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
+    OVector<T, D>: Clone,
+    OMatrix<T, D, D>: Clone,
+{
+    fn clone(&self) -> Self {
+        Self {
+            eigenvalues: self.eigenvalues.clone(),
+            eigenvectors: self.eigenvectors.clone(),
+            left_eigenvectors: self.left_eigenvectors.clone(),
+        }
+    }
+}
+
+impl<T: fmt::Debug, D: Dim> fmt::Debug for Eigen<T, D>
+where
+    DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
+    OVector<T, D>: fmt::Debug,
+    OMatrix<T, D, D>: fmt::Debug,
+{
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        f.debug_struct("Eigen")
+            .field("eigenvalues", &self.eigenvalues)
+            .field("eigenvectors", &self.eigenvectors)
+            .field("left_eigenvectors", &self.left_eigenvectors)
+            .finish()
+    }
+}
+
 impl<T: EigenScalar + RealField, D: Dim> Eigen<T, D>
 where
     DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
diff --git a/src/base/allocator.rs b/src/base/allocator.rs
index 95a65c6f..26ea11bc 100644
--- a/src/base/allocator.rs
+++ b/src/base/allocator.rs
@@ -17,7 +17,8 @@ use crate::base::DefaultAllocator;
 /// Every allocator must be both static and dynamic. Though not all implementations may share the
 /// same `Buffer` type.
 ///
-/// If you also want to be able to create uninitizalized memory buffers, see [`Allocator`].
+/// If you also want to be able to create uninitizalized or manually dropped memory buffers, see
+/// [`Allocator`].
 pub trait InnerAllocator<T, R: Dim, C: Dim = U1>: 'static + Sized {
     /// The type of buffer this allocator can instanciate.
     type Buffer: ContiguousStorageMut<T, R, C>;
@@ -44,6 +45,10 @@ pub trait Allocator<T, R: Dim, C: Dim = U1>:
     ) -> <Self as InnerAllocator<MaybeUninit<T>, R, C>>::Buffer;
 
     /// Assumes a data buffer to be initialized. This operation should be near zero-cost.
+    ///
+    /// # Safety
+    /// The user must make sure that every single entry of the buffer has been initialized,
+    /// or Undefined Behavior will immediately occur.    
     unsafe fn assume_init(
         uninit: <Self as InnerAllocator<MaybeUninit<T>, R, C>>::Buffer,
     ) -> <Self as InnerAllocator<T, R, C>>::Buffer;
diff --git a/src/base/array_storage.rs b/src/base/array_storage.rs
index b87442a4..bcf9df33 100644
--- a/src/base/array_storage.rs
+++ b/src/base/array_storage.rs
@@ -1,4 +1,4 @@
-use std::fmt::{self, Debug, Formatter};
+use std::mem;use std::fmt::{self, Debug, Formatter};
 // use std::hash::{Hash, Hasher};
 #[cfg(feature = "abomonation-serialize")]
 use std::io::{Result as IOResult, Write};
@@ -31,7 +31,7 @@ use crate::base::storage::{
  *
  */
 /// A array-based statically sized matrix data storage.
-#[repr(C)]
+#[repr(transparent)]
 #[derive(Copy, Clone, PartialEq, Eq, Hash)]
 pub struct ArrayStorage<T, const R: usize, const C: usize>(pub [[T; R]; C]);
 
@@ -155,8 +155,8 @@ where
 
     fn reshape_generic(self, _: Const<R2>, _: Const<C2>) -> Self::Output {
         unsafe {
-            let data: [[T; R2]; C2] = std::mem::transmute_copy(&self.0);
-            std::mem::forget(self.0);
+            let data: [[T; R2]; C2] = mem::transmute_copy(&self.0);
+            mem::forget(self.0);
             ArrayStorage(data)
         }
     }
diff --git a/src/base/blas.rs b/src/base/blas.rs
index 9654df08..4f605e0f 100644
--- a/src/base/blas.rs
+++ b/src/base/blas.rs
@@ -6,7 +6,7 @@
 //! that return an owned matrix that would otherwise result from setting a
 //! parameter to zero in the other methods.
 
-use crate::SimdComplexField;
+use crate::{MatrixSliceMut, SimdComplexField, VectorSliceMut};
 #[cfg(feature = "std")]
 use matrixmultiply;
 use num::{One, Zero};
@@ -717,10 +717,15 @@ where
     /// Computes `alpha * a * x`, where `a` is a matrix, `x` a vector, and
     /// `alpha` is a scalar.
     ///
-    /// # Safety
     /// `self` must be completely uninitialized, or data leaks will occur. After
     /// this method is called, all entries in `self` will be initialized.
-    pub fn axc<D2: Dim, S2>(&mut self, a: T, x: &Vector<T, D2, S2>, c: T)
+    #[inline]
+    pub fn axc<D2: Dim, S2>(
+        &mut self,
+        a: T,
+        x: &Vector<T, D2, S2>,
+        c: T,
+    ) -> VectorSliceMut<T, D, S::RStride, S::CStride>
     where
         S2: Storage<T, D2>,
         ShapeConstraint: DimEq<D, D2>,
@@ -728,10 +733,15 @@ where
         let rstride1 = self.strides().0;
         let rstride2 = x.strides().0;
 
+        // Safety: see each individual remark.
         unsafe {
+            // We don't mind `x` and `y` not being contiguous, as we'll only
+            // access the elements we're allowed to. (TODO: double check this)
             let y = self.data.as_mut_slice_unchecked();
             let x = x.data.as_slice_unchecked();
 
+            // The indices are within range, and only access elements that belong
+            // to `x` and `y` themselves.
             for i in 0..y.len() {
                 *y.get_unchecked_mut(i * rstride1) = MaybeUninit::new(
                     a.inlined_clone()
@@ -739,20 +749,26 @@ where
                         * c.inlined_clone(),
                 );
             }
+
+            // We've initialized all elements.
+            self.assume_init_mut()
         }
     }
 
     /// Computes `alpha * a * x`, where `a` is a matrix, `x` a vector, and
     /// `alpha` is a scalar.
     ///
-    /// Initializes `self`.
+    /// `self` must be completely uninitialized, or data leaks will occur. After
+    /// the method is called, `self` will be completely initialized. We return
+    /// an initialized mutable vector slice to `self` for convenience.
     #[inline]
     pub fn gemv_z<R2: Dim, C2: Dim, D3: Dim, SB, SC>(
         &mut self,
         alpha: T,
         a: &Matrix<T, R2, C2, SB>,
         x: &Vector<T, D3, SC>,
-    ) where
+    ) -> VectorSliceMut<T, D, S::RStride, S::CStride>
+    where
         T: One,
         SB: Storage<T, R2, C2>,
         SC: Storage<T, D3>,
@@ -769,24 +785,28 @@ where
 
         if ncols2 == 0 {
             self.fill_fn(|| MaybeUninit::new(T::zero()));
-            return;
+
+            // Safety: all entries have just been initialized.
+            unsafe {
+                return self.assume_init_mut();
+            }
         }
 
         // 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 init = self.axc(alpha.inlined_clone(), &col2, val);
 
-        // Safety: axc initializes self.
+        // Safety: all indices are within range.
         unsafe {
-            let mut init = self.assume_init_mut();
-
             for j in 1..ncols2 {
                 let col2 = a.column(j);
                 let val = x.vget_unchecked(j).inlined_clone();
                 init.axcpy(alpha.inlined_clone(), &col2, val, T::one());
             }
         }
+
+        init
     }
 
     #[inline(always)]
@@ -825,9 +845,8 @@ where
         // 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 = 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 {
@@ -894,7 +913,8 @@ where
         alpha: T,
         a: &Matrix<T, R2, C2, SB>,
         b: &Matrix<T, R3, C3, SC>,
-    ) where
+    ) -> MatrixSliceMut<T, R1, C1, S::RStride, S::CStride>
+    where
         SB: Storage<T, R2, C2>,
         SC: Storage<T, R3, C3>,
         ShapeConstraint: SameNumberOfRows<R1, R2>
@@ -945,7 +965,9 @@ where
                     // enter this codepath.
                     if ncols1 == 0 {
                         self.fill_fn(|| MaybeUninit::new(T::zero()));
-                        return;
+
+                        // Safety: there's no (uninitialized) values.
+                        return unsafe{self.assume_init_mut()};
                     }
 
                     let (rsa, csa) = a.strides();
@@ -970,8 +992,6 @@ where
                                 rsc as isize,
                                 csc as isize,
                             );
-
-                            return;
                         }
                     } else if T::is::<f64>() {
                         unsafe {
@@ -991,19 +1011,26 @@ where
                                 rsc as isize,
                                 csc as isize,
                             );
-
-                            return;
                         }
                     }
+
+                    // Safety: all entries have been initialized.
+                    unsafe {
+                        return self.assume_init_mut();
+                    }
                 }
             }
         }
 
         for j1 in 0..ncols1 {
             // TODO: avoid bound checks.
-            self.column_mut(j1)
+            let _ = self
+                .column_mut(j1)
                 .gemv_z(alpha.inlined_clone(), a, &b.column(j1));
         }
+
+        // Safety: all entries have been initialized.
+        unsafe { self.assume_init_mut() }
     }
 }
 
@@ -1571,8 +1598,7 @@ where
     {
         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 mut work = unsafe { work.assume_init() };
+        let mut work = work.gemv_z(T::one(), lhs, &mid.column(0));
 
         self.ger(alpha.inlined_clone(), &work, &lhs.column(0), beta);
 
@@ -1614,14 +1640,12 @@ where
     ) where
         S3: Storage<T, D3, D3>,
         S4: Storage<T, R4, C4>,
-        ShapeConstraint: DimEq<R4, D3> + DimEq<D3, R4> + DimEq<D1, C4>,
+        ShapeConstraint: DimEq<D3, R4> + DimEq<R4, D3> + DimEq<D1, C4>,
         DefaultAllocator: Allocator<T, D3>,
     {
         // TODO: figure out why type inference isn't doing its job.
-        let mut work =
-            Matrix::new_uninitialized_generic(D3::from_usize(mid.shape().0), Const::<1>);
-        work.gemv_z::<D3, D3, R4, S3, _>(T::one(), mid, &rhs.column(0));
-        let mut work = unsafe { work.assume_init() };
+        let mut work = Matrix::new_uninitialized_generic(D3::from_usize(mid.shape().0), Const::<1>);
+        let mut work = work.gemv_z::<D3, _, _, _, _>(T::one(), mid, &rhs.column(0));
 
         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 b768ed73..b8a50048 100644
--- a/src/base/conversion.rs
+++ b/src/base/conversion.rs
@@ -1,9 +1,10 @@
+use std::borrow::{Borrow, BorrowMut};
+use std::convert::{AsMut, AsRef, From, Into};
+use std::mem::{self, ManuallyDrop, MaybeUninit};
+
 #[cfg(all(feature = "alloc", not(feature = "std")))]
 use alloc::vec::Vec;
 use simba::scalar::{SubsetOf, SupersetOf};
-use std::borrow::{Borrow, BorrowMut};
-use std::convert::{AsMut, AsRef, From, Into};
-use std::mem::MaybeUninit;
 
 use simba::simd::{PrimitiveSimdValue, SimdValue};
 
@@ -105,18 +106,18 @@ impl<'a, T, R: Dim, C: Dim, S: StorageMut<T, R, C>> IntoIterator for &'a mut Mat
 impl<T, const D: usize> From<[T; D]> for SVector<T, D> {
     #[inline]
     fn from(arr: [T; D]) -> Self {
-        unsafe { Self::from_data_statically_unchecked(ArrayStorage([arr; 1])) }
+        Self::from_data(ArrayStorage([arr; 1]))
     }
 }
 
-impl<T, const D: usize> From<SVector<T, D>> for [T; D]
-where
-    T: Clone,
-{
+impl<T, const D: usize> From<SVector<T, D>> for [T; D] {
     #[inline]
     fn from(vec: SVector<T, D>) -> Self {
-        // TODO: unfortunately, we must clone because we can move out of an array.
-        vec.data.0[0].clone()
+        let data = ManuallyDrop::new(vec.data.0);
+        // Safety: [[T; D]; 1] always has the same data layout as [T; D].
+        let res = unsafe { (data.as_ptr() as *const [_; D]).read() };
+        mem::forget(data);
+        res
     }
 }
 
@@ -184,7 +185,7 @@ impl_from_into_asref_1D!(
 impl<T, const R: usize, const C: usize> From<[[T; R]; C]> for SMatrix<T, R, C> {
     #[inline]
     fn from(arr: [[T; R]; C]) -> Self {
-        unsafe { Self::from_data_statically_unchecked(ArrayStorage(arr)) }
+        Self::from_data(ArrayStorage(arr))
     }
 }
 
@@ -326,7 +327,8 @@ where
                 (row_slice, col_slice),
                 (rstride_slice, cstride_slice),
             );
-            Matrix::from_data_statically_unchecked(data)
+
+            Self::from_data(data)
         }
     }
 }
@@ -356,7 +358,8 @@ where
                 (row_slice, col_slice),
                 (rstride_slice, cstride_slice),
             );
-            Matrix::from_data_statically_unchecked(data)
+
+            Matrix::from_data(data)
         }
     }
 }
@@ -386,7 +389,8 @@ where
                 (row_slice, col_slice),
                 (rstride_slice, cstride_slice),
             );
-            Matrix::from_data_statically_unchecked(data)
+
+            Matrix::from_data(data)
         }
     }
 }
diff --git a/src/base/default_allocator.rs b/src/base/default_allocator.rs
index 4d8d0010..b30e8960 100644
--- a/src/base/default_allocator.rs
+++ b/src/base/default_allocator.rs
@@ -76,11 +76,10 @@ impl<T, const R: usize, const C: usize> Allocator<T, Const<R>, Const<C>> for Def
     unsafe fn assume_init(
         uninit: <Self as InnerAllocator<MaybeUninit<T>, Const<R>, Const<C>>>::Buffer,
     ) -> Owned<T, Const<R>, Const<C>> {
-        // SAFETY:
+        // Safety:
         // * The caller guarantees that all elements of the array are initialized
         // * `MaybeUninit<T>` and T are guaranteed to have the same layout
         // * `MaybeUnint` does not drop, so there are no double-frees
-        // * `ArrayStorage` is transparent.
         // And thus the conversion is safe
         ArrayStorage((&uninit as *const _ as *const [_; C]).read())
     }
diff --git a/src/base/dimension.rs b/src/base/dimension.rs
index 8573dd59..22b80b2a 100644
--- a/src/base/dimension.rs
+++ b/src/base/dimension.rs
@@ -2,7 +2,7 @@
 
 //! Traits and tags for identifying the dimension of all algebraic entities.
 
-use std::any::{Any, TypeId};
+use std::any::TypeId;
 use std::cmp;
 use std::fmt::Debug;
 use std::ops::{Add, Div, Mul, Sub};
@@ -11,7 +11,7 @@ use typenum::{self, Diff, Max, Maximum, Min, Minimum, Prod, Quot, Sum, Unsigned}
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
 
-/// Dim of dynamically-sized algebraic entities.
+/// Stores the dimension of dynamically-sized algebraic entities.
 #[derive(Clone, Copy, Eq, PartialEq, Debug)]
 pub struct Dynamic {
     value: usize,
@@ -55,7 +55,7 @@ impl IsNotStaticOne for Dynamic {}
 
 /// Trait implemented by any type that can be used as a dimension. This includes type-level
 /// integers and `Dynamic` (for dimensions not known at compile-time).
-pub trait Dim: Any + Debug + Copy + PartialEq + Send + Sync {
+pub trait Dim: 'static + Debug + Copy + PartialEq + Send + Sync {
     #[inline(always)]
     fn is<D: Dim>() -> bool {
         TypeId::of::<Self>() == TypeId::of::<D>()
@@ -196,6 +196,9 @@ dim_ops!(
     DimMax, DimNameMax, Max, max, cmp::max, DimMaximum, DimNameMaximum, Maximum;
 );
 
+/// A wrapper around const types, which provides the capability of performing
+/// type-level arithmetic. This might get removed if const-generics become
+/// more powerful in the future.
 #[derive(Debug, Copy, Clone, PartialEq, Eq, Hash)]
 pub struct Const<const R: usize>;
 
diff --git a/src/base/indexing.rs b/src/base/indexing.rs
index 0073c85f..a8db21ec 100644
--- a/src/base/indexing.rs
+++ b/src/base/indexing.rs
@@ -673,7 +673,7 @@ macro_rules! impl_index_pair {
                         (rows.lower(nrows),  cols.lower(ncols)),
                         (rows.length(nrows), cols.length(ncols)));
 
-                Matrix::from_data_statically_unchecked(data)
+                Matrix::from_data(data)
             }
         }
 
@@ -699,7 +699,7 @@ macro_rules! impl_index_pair {
                         (rows.lower(nrows),  cols.lower(ncols)),
                         (rows.length(nrows), cols.length(ncols)));
 
-                Matrix::from_data_statically_unchecked(data)
+                Matrix::from_data(data)
             }
         }
     }
diff --git a/src/base/matrix.rs b/src/base/matrix.rs
index 6ef2c162..94c3f88e 100644
--- a/src/base/matrix.rs
+++ b/src/base/matrix.rs
@@ -5,10 +5,11 @@ use std::io::{Result as IOResult, Write};
 use approx::{AbsDiffEq, RelativeEq, UlpsEq};
 use std::any::TypeId;
 use std::cmp::Ordering;
-use std::fmt;use std::ptr;
+use std::fmt;
 use std::hash::{Hash, Hasher};
 use std::marker::PhantomData;
 use std::mem::{self, ManuallyDrop, MaybeUninit};
+use std::ptr;
 
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
@@ -26,7 +27,7 @@ use crate::base::iter::{
     ColumnIter, ColumnIterMut, MatrixIter, MatrixIterMut, RowIter, RowIterMut,
 };
 use crate::base::storage::{
-    ContiguousStorage, ContiguousStorageMut, Owned, SameShapeStorage, Storage, StorageMut,
+    ContiguousStorage, ContiguousStorageMut, SameShapeStorage, Storage, StorageMut,
 };
 use crate::base::{Const, DefaultAllocator, OMatrix, OVector, Scalar, Unit};
 use crate::{ArrayStorage, MatrixSlice, MatrixSliceMut, SMatrix, SimdComplexField};
@@ -151,7 +152,7 @@ pub type MatrixCross<T, R1, C1, R2, C2> =
 /// Note that mixing `Dynamic` with type-level unsigned integers is allowed. Actually, a
 /// 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)]
+#[repr(transparent)]
 #[derive(Clone, Copy)]
 pub struct Matrix<T, R, C, S> {
     /// The data storage that contains all the matrix components. Disappointed?
@@ -187,8 +188,8 @@ pub struct Matrix<T, R, C, S> {
     //       Note that it would probably make sense to just have
     //       the type `Matrix<S>`, and have `T, R, C` be associated-types
     //       of the `Storage` trait. However, because we don't have
-    //       specialization, this is not bossible because these `T, R, C`
-    //       allows us to desambiguate a lot of configurations.
+    //       specialization, this is not possible because these `T, R, C`
+    //       allows us to disambiguate a lot of configurations.
     _phantoms: PhantomData<(T, R, C)>,
 }
 
@@ -198,9 +199,12 @@ impl<T, R: Dim, C: Dim, S: fmt::Debug> fmt::Debug for Matrix<T, R, C, S> {
     }
 }
 
-impl<T, R: Dim, C: Dim, S: Default> Default for Matrix<T, R, C, S> {
+impl<T, R: Dim, C: Dim, S> Default for Matrix<T, R, C, S>
+where
+    S: Storage<T, R, C> + Default,
+{
     fn default() -> Self {
-        unsafe { Matrix::from_data_statically_unchecked(Default::default()) }
+        Matrix::from_data(Default::default())
     }
 }
 
@@ -330,8 +334,19 @@ mod rkyv_impl {
 }
 
 impl<T, R, C, S> Matrix<T, R, C, S> {
-    /// Creates a new matrix with the given data without statically checking that the matrix
-    /// dimension matches the storage dimension.
+    /// Creates a new matrix with the given data without statically checking
+    /// that the matrix dimension matches the storage dimension.
+    ///
+    /// There's only two instances in which you should use this method instead
+    /// of the safe counterpart [`from_data`]:
+    /// - You can't get the type checker to validate your matrices, even though
+    ///   you're **certain** that they're of the right dimensions.
+    /// - You want to declare a matrix in a `const` context.
+    ///
+    /// # Safety
+    /// If the storage dimension does not match the matrix dimension, any other
+    /// method called on this matrix may behave erroneously, panic, or cause
+    /// Undefined Behavior.
     #[inline(always)]
     pub const unsafe fn from_data_statically_unchecked(data: S) -> Matrix<T, R, C, S> {
         Matrix {
@@ -348,21 +363,17 @@ where
 {
     /// Allocates a matrix with the given number of rows and columns without initializing its content.
     pub fn new_uninitialized_generic(nrows: R, ncols: C) -> OMatrix<MaybeUninit<T>, R, C> {
-        unsafe {
-            OMatrix::from_data_statically_unchecked(
-                <DefaultAllocator as Allocator<T, R, C>>::allocate_uninitialized(nrows, ncols),
-            )
-        }
+        OMatrix::from_data(
+            <DefaultAllocator as Allocator<T, R, C>>::allocate_uninitialized(nrows, ncols),
+        )
     }
 
     /// Converts this matrix into one whose entries need to be manually dropped. This should be
     /// near zero-cost.
     pub fn manually_drop(self) -> OMatrix<ManuallyDrop<T>, R, C> {
-        unsafe {
-            OMatrix::from_data_statically_unchecked(
-                <DefaultAllocator as Allocator<T, R, C>>::manually_drop(self.data),
-            )
-        }
+        OMatrix::from_data(<DefaultAllocator as Allocator<T, R, C>>::manually_drop(
+            self.data,
+        ))
     }
 }
 
@@ -375,19 +386,21 @@ where
     ///
     /// For the similar method that operates on matrix slices, see [`slice_assume_init`].
     pub unsafe fn assume_init(self) -> OMatrix<T, R, C> {
-        OMatrix::from_data_statically_unchecked(
-            <DefaultAllocator as Allocator<T, R, C>>::assume_init(self.data),
-        )
+        OMatrix::from_data(<DefaultAllocator as Allocator<T, R, C>>::assume_init(
+            self.data,
+        ))
     }
 
-    /// Assumes a matrix's entries to be initialized, and drops them. This allows the
-    /// buffer to be safely reused.
-    pub fn reinitialize(&mut self) {
+    /// Assumes a matrix's entries to be initialized, and drops them in place.
+    /// This allows the buffer to be safely reused.
+    ///
+    /// # Safety
+    /// All of the matrix's entries need to be uninitialized. Otherwise,
+    /// Undefined Behavior will be triggered.
+    pub unsafe fn reinitialize(&mut self) {
         for i in 0..self.nrows() {
             for j in 0..self.ncols() {
-                unsafe {
-                    ptr::drop_in_place(self.get_unchecked_mut((i, j)));
-                }
+                ptr::drop_in_place(self.get_unchecked_mut((i, j)));
             }
         }
     }
@@ -418,8 +431,8 @@ impl<T, const R: usize, const C: usize> SMatrix<T, R, C> {
     /// work in `const fn` contexts.
     #[inline(always)]
     pub const fn from_array_storage(storage: ArrayStorage<T, R, C>) -> Self {
-        // This is sound because the row and column types are exactly the same as that of the
-        // storage, so there can be no mismatch
+        // Safety: This is sound because the row and column types are exactly
+        // the same as that of the storage, so there can be no mismatch.
         unsafe { Self::from_data_statically_unchecked(storage) }
     }
 }
@@ -433,8 +446,8 @@ impl<T> DMatrix<T> {
     /// This method exists primarily as a workaround for the fact that `from_data` can not
     /// work in `const fn` contexts.
     pub const fn from_vec_storage(storage: VecStorage<T, Dynamic, Dynamic>) -> Self {
-        // This is sound because the dimensions of the matrix and the storage are guaranteed
-        // to be the same
+        // Safety: This is sound because the dimensions of the matrix and the
+        // storage are guaranteed to be the same.
         unsafe { Self::from_data_statically_unchecked(storage) }
     }
 }
@@ -448,8 +461,8 @@ impl<T> DVector<T> {
     /// This method exists primarily as a workaround for the fact that `from_data` can not
     /// work in `const fn` contexts.
     pub const fn from_vec_storage(storage: VecStorage<T, Dynamic, U1>) -> Self {
-        // This is sound because the dimensions of the matrix and the storage are guaranteed
-        // to be the same
+        // Safety: This is sound because the dimensions of the matrix and the
+        // storage are guaranteed to be the same.
         unsafe { Self::from_data_statically_unchecked(storage) }
     }
 }
@@ -458,6 +471,8 @@ impl<T, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Creates a new matrix with the given data.
     #[inline(always)]
     pub fn from_data(data: S) -> Self {
+        // Safety: This is sound because the dimensions of the matrix and the
+        // storage are guaranteed to be the same.
         unsafe { Self::from_data_statically_unchecked(data) }
     }
 
@@ -623,19 +638,22 @@ impl<T, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     #[inline]
     pub fn into_owned_sum<R2: Dim, C2: Dim>(self) -> MatrixSum<T, R, C, R2, C2>
     where
-        T: Clone + 'static,
+        T: Clone,
         DefaultAllocator: SameShapeAllocator<T, R, C, R2, C2>,
         ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
     {
-        if TypeId::of::<SameShapeStorage<T, R, C, R2, C2>>() == TypeId::of::<Owned<T, R, C>>() {
-            // We can just return `self.into_owned()`.
-
+        // If both storages are the same, we can just return `self.into_owned()`.
+        // Unfortunately, it's not trivial to convince the compiler of this.
+        if TypeId::of::<SameShapeR<R, R2>>() == TypeId::of::<R>()
+            && TypeId::of::<SameShapeC<C, C2>>() == TypeId::of::<C>()
+        {
+            // Safety: we're transmuting from a type into itself, and we make
+            // sure not to leak anything.
             unsafe {
-                // TODO: check that those copies are optimized away by the compiler.
-                let owned = self.into_owned();
-                let res = mem::transmute_copy(&owned);
-                mem::forget(owned);
-                res
+                let mat = self.into_owned();
+                let mat_copy = mem::transmute_copy(&mat);
+                mem::forget(mat);
+                mat_copy
             }
         } else {
             self.clone_owned_sum()
diff --git a/src/base/matrix_slice.rs b/src/base/matrix_slice.rs
index 7ba2eb8d..25baee55 100644
--- a/src/base/matrix_slice.rs
+++ b/src/base/matrix_slice.rs
@@ -222,7 +222,12 @@ storage_impl!(SliceStorage, SliceStorageMut);
 impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim>
     SliceStorage<'a, MaybeUninit<T>, R, C, RStride, CStride>
 {
-    /// Assumes a slice storage's entries to be initialized. This operation should be near zero-cost.
+    /// Assumes a slice storage's entries to be initialized. This operation
+    /// should be near zero-cost.
+    ///
+    /// # Safety
+    /// All of the slice storage's entries must be initialized, otherwise
+    /// Undefined Behavior will be triggered.
     pub unsafe fn assume_init(self) -> SliceStorage<'a, T, R, C, RStride, CStride> {
         SliceStorage::from_raw_parts(self.ptr as *const T, self.shape, self.strides)
     }
@@ -401,7 +406,7 @@ macro_rules! matrix_slice_impl(
 
             unsafe {
                 let data = $SliceStorage::new_unchecked($data, (row_start, 0), shape);
-                Matrix::from_data_statically_unchecked(data)
+                Matrix::from_data(data)
             }
         }
 
@@ -421,7 +426,7 @@ macro_rules! matrix_slice_impl(
 
             unsafe {
                 let data = $SliceStorage::new_with_strides_unchecked($data, (row_start, 0), shape, strides);
-                Matrix::from_data_statically_unchecked(data)
+                Matrix::from_data(data)
             }
         }
 
@@ -488,7 +493,7 @@ macro_rules! matrix_slice_impl(
 
             unsafe {
                 let data = $SliceStorage::new_unchecked($data, (0, first_col), shape);
-                Matrix::from_data_statically_unchecked(data)
+                Matrix::from_data(data)
             }
         }
 
@@ -508,7 +513,7 @@ macro_rules! matrix_slice_impl(
 
             unsafe {
                 let data = $SliceStorage::new_with_strides_unchecked($data, (0, first_col), shape, strides);
-                Matrix::from_data_statically_unchecked(data)
+                Matrix::from_data(data)
             }
         }
 
@@ -528,7 +533,7 @@ macro_rules! matrix_slice_impl(
 
             unsafe {
                 let data = $SliceStorage::new_unchecked($data, start, shape);
-                Matrix::from_data_statically_unchecked(data)
+                Matrix::from_data(data)
             }
         }
 
@@ -555,7 +560,7 @@ macro_rules! matrix_slice_impl(
 
             unsafe {
                 let data = $SliceStorage::new_unchecked($data, (irow, icol), shape);
-                Matrix::from_data_statically_unchecked(data)
+                Matrix::from_data(data)
             }
         }
 
@@ -579,7 +584,7 @@ macro_rules! matrix_slice_impl(
 
             unsafe {
                 let data = $SliceStorage::new_unchecked($data, start, shape);
-                Matrix::from_data_statically_unchecked(data)
+                Matrix::from_data(data)
             }
         }
 
@@ -601,7 +606,7 @@ macro_rules! matrix_slice_impl(
 
             unsafe {
                 let data = $SliceStorage::new_with_strides_unchecked($data, start, shape, strides);
-                Matrix::from_data_statically_unchecked(data)
+                Matrix::from_data(data)
             }
         }
 
@@ -645,8 +650,8 @@ macro_rules! matrix_slice_impl(
 
                 let data1  = $SliceStorage::from_raw_parts(ptr1, (nrows1, ncols), strides);
                 let data2  = $SliceStorage::from_raw_parts(ptr2, (nrows2, ncols), strides);
-                let slice1 = Matrix::from_data_statically_unchecked(data1);
-                let slice2 = Matrix::from_data_statically_unchecked(data2);
+                let slice1 = Matrix::from_data(data1);
+                let slice2 = Matrix::from_data(data2);
 
                 (slice1, slice2)
             }
@@ -681,8 +686,8 @@ macro_rules! matrix_slice_impl(
 
                 let data1  = $SliceStorage::from_raw_parts(ptr1, (nrows, ncols1), strides);
                 let data2  = $SliceStorage::from_raw_parts(ptr2, (nrows, ncols2), strides);
-                let slice1 = Matrix::from_data_statically_unchecked(data1);
-                let slice2 = Matrix::from_data_statically_unchecked(data2);
+                let slice1 = Matrix::from_data(data1);
+                let slice2 = Matrix::from_data(data2);
 
                 (slice1, slice2)
             }
@@ -1007,6 +1012,6 @@ impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim>
             _phantoms: PhantomData,
         };
 
-        unsafe { Matrix::from_data_statically_unchecked(data) }
+       Matrix::from_data(data) 
     }
 }
diff --git a/src/base/ops.rs b/src/base/ops.rs
index 25921e90..dfedb69a 100644
--- a/src/base/ops.rs
+++ b/src/base/ops.rs
@@ -17,7 +17,7 @@ use crate::base::dimension::{Dim, DimMul, DimName, DimProd, Dynamic};
 use crate::base::storage::{ContiguousStorageMut, Storage, StorageMut};
 use crate::base::{DefaultAllocator, Matrix, MatrixSum, OMatrix, Scalar, VectorSlice};
 use crate::storage::Owned;
-use crate::SimdComplexField;
+use crate::{MatrixSliceMut, SimdComplexField};
 
 /*
  *
@@ -581,7 +581,7 @@ where
     #[inline]
     fn mul(self, rhs: &'b Matrix<T, R2, C2, SB>) -> Self::Output {
         let mut res = Matrix::new_uninitialized_generic(self.data.shape().0, rhs.data.shape().1);
-        self.mul_to(rhs, &mut res);
+        let _ = self.mul_to(rhs, &mut res);
         unsafe { res.assume_init() }
     }
 }
@@ -645,7 +645,7 @@ impl<T, R1: Dim, C1: Dim, R2: Dim, SA, SB> MulAssign<Matrix<T, R2, C1, SB>>
 where
     T: Scalar + Zero + One + ClosedAdd + ClosedMul,
     SB: Storage<T, R2, C1>,
-    SA: ContiguousStorageMut<T, R1, C1> ,
+    SA: ContiguousStorageMut<T, R1, C1>,
     ShapeConstraint: AreMultipliable<R1, C1, R2, C1>,
     DefaultAllocator: Allocator<T, R1, C1> + InnerAllocator<T, R1, C1, Buffer = SA>,
 {
@@ -660,7 +660,7 @@ impl<'b, T, R1: Dim, C1: Dim, R2: Dim, SA, SB> MulAssign<&'b Matrix<T, R2, C1, S
 where
     T: Scalar + Zero + One + ClosedAdd + ClosedMul,
     SB: Storage<T, R2, C1>,
-    SA: ContiguousStorageMut<T, R1, C1> ,
+    SA: ContiguousStorageMut<T, R1, C1>,
     ShapeConstraint: AreMultipliable<R1, C1, R2, C1>,
     // TODO: this is too restrictive. See comments for the non-ref version.
     DefaultAllocator: Allocator<T, R1, C1> + InnerAllocator<T, R1, C1, Buffer = SA>,
@@ -786,18 +786,19 @@ where
 
     /// Equivalent to `self * rhs` but stores the result into `out` to avoid allocations.
     #[inline]
-    pub fn mul_to<R2: Dim, C2: Dim, SB, R3: Dim, C3: Dim, SC>(
+    pub fn mul_to<'a, R2: Dim, C2: Dim, SB, R3: Dim, C3: Dim, SC>(
         &self,
         rhs: &Matrix<T, R2, C2, SB>,
-        out: &mut Matrix<MaybeUninit<T>, R3, C3, SC>,
-    ) where
+        out: &'a mut Matrix<MaybeUninit<T>, R3, C3, SC>,
+    ) -> MatrixSliceMut<'a, T, R3, C3, SC::RStride, SC::CStride>
+    where
         SB: Storage<T, R2, C2>,
         SC: StorageMut<MaybeUninit<T>, R3, C3>,
         ShapeConstraint: SameNumberOfRows<R3, R1>
             + SameNumberOfColumns<C3, C2>
             + AreMultipliable<R1, C1, R2, C2>,
     {
-        out.gemm_z(T::one(), self, rhs);
+        out.gemm_z(T::one(), self, rhs)
     }
 
     /// The kronecker product of two matrices (aka. tensor product of the corresponding linear
diff --git a/src/base/scalar.rs b/src/base/scalar.rs
index c14f3eb7..80a78594 100644
--- a/src/base/scalar.rs
+++ b/src/base/scalar.rs
@@ -10,20 +10,24 @@ use std::fmt::Debug;
 /// - Makes debugging generic code possible in most circumstances.
 pub trait Scalar: 'static + Clone + Debug {
     #[inline]
-    /// Tests if `Self` is the same as the type `T`.
+    /// Tests whether `Self` is the same as the type `T`.
     ///
     /// Typically used to test of `Self` is an `f32` or an `f64`, which is
     /// important as it allows for specialization and certain optimizations to
     /// be made.
     ///
-    /// If the need ever arose to get rid of the `'static` requirement
+    // If the need ever arose to get rid of the `'static` requirement, we could
+    // merely replace this method by two unsafe associated methods `is_f32` and
+    // `is_f64`.
     fn is<T: Scalar>() -> bool {
         TypeId::of::<Self>() == TypeId::of::<T>()
     }
 
     /// Performance hack: Clone doesn't get inlined for Copy types in debug
     /// mode, so make it inline anyway.
-    fn inlined_clone(&self) -> Self;
+    fn inlined_clone(&self) -> Self {
+        self.clone()
+    }
 }
 
 // Unfortunately, this blanket impl leads to many misleading compiler messages
diff --git a/src/base/unit.rs b/src/base/unit.rs
index f656b247..ed9ffc14 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 _ as *const Self) }
+        unsafe { &*(value as *const _ as *const _) }
     }
 
     /// Retrieves the underlying value.
diff --git a/src/base/vec_storage.rs b/src/base/vec_storage.rs
index ee57218f..9f9d649d 100644
--- a/src/base/vec_storage.rs
+++ b/src/base/vec_storage.rs
@@ -28,7 +28,6 @@ use abomonation::Abomonation;
  *
  */
 /// A Vec-based matrix data storage. It may be dynamically-sized.
-#[repr(C)]
 #[derive(Eq, Debug, Clone, PartialEq)]
 pub struct VecStorage<T, R: Dim, C: Dim> {
     pub(crate) data: Vec<T>,
diff --git a/src/geometry/dual_quaternion.rs b/src/geometry/dual_quaternion.rs
index 17af51fe..0469829f 100644
--- a/src/geometry/dual_quaternion.rs
+++ b/src/geometry/dual_quaternion.rs
@@ -279,6 +279,7 @@ impl<'a, T: Deserialize<'a>> Deserialize<'a> for DualQuaternion<T> {
 
 impl<T> DualQuaternion<T> {
     // TODO: Cloning shouldn't be necessary.
+    // TODO: rename into `into_vector` to appease clippy.
     fn to_vector(self) -> OVector<T, U8>
     where
         T: Clone,
diff --git a/src/geometry/dual_quaternion_ops.rs b/src/geometry/dual_quaternion_ops.rs
index 4f1e58e3..151b2e05 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 _ as *const [T; 8]) }
+        unsafe { &*(self as *const _ as *const _) }
     }
 }
 
 impl<T> AsMut<[T; 8]> for DualQuaternion<T> {
     #[inline]
     fn as_mut(&mut self) -> &mut [T; 8] {
-        unsafe { &mut *(self as *mut _ as *mut [T; 8]) }
+        unsafe { &mut *(self as *mut _ as *mut _) }
     }
 }
 
diff --git a/src/geometry/isometry.rs b/src/geometry/isometry.rs
index cb56ad83..389965be 100755
--- a/src/geometry/isometry.rs
+++ b/src/geometry/isometry.rs
@@ -53,7 +53,6 @@ use crate::geometry::{AbstractRotation, Point, Translation};
 /// # Conversion to a matrix
 /// * [Conversion to a matrix <span style="float:right;">`to_matrix`…</span>](#conversion-to-a-matrix)
 ///
-#[repr(C)]
 #[derive(Debug)]
 #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))]
 #[cfg_attr(
diff --git a/src/geometry/orthographic.rs b/src/geometry/orthographic.rs
index 974df3ff..ba613de7 100644
--- a/src/geometry/orthographic.rs
+++ b/src/geometry/orthographic.rs
@@ -18,7 +18,7 @@ use crate::base::{Matrix4, Vector, Vector3};
 use crate::geometry::{Point3, Projective3};
 
 /// A 3D orthographic projection stored as a homogeneous 4x4 matrix.
-#[repr(C)]
+#[repr(transparent)]
 pub struct Orthographic3<T> {
     matrix: Matrix4<T>,
 }
@@ -235,6 +235,7 @@ impl<T> Orthographic3<T> {
     /// ```
     #[inline]
     #[must_use]
+    // TODO: rename into `into_homogeneous` to appease clippy.
     pub fn to_homogeneous(self) -> Matrix4<T> {
         self.matrix
     }
@@ -270,8 +271,8 @@ impl<T> Orthographic3<T> {
     #[inline]
     #[must_use]
     pub fn as_projective(&self) -> &Projective3<T> {
-        // Safety: Self and Projective3 are both #[repr(C)] of a matrix.
-        unsafe { &*(self as *const _ as *const Projective3<T>) }
+        // Safety: Self and Projective3 are both #[repr(transparent)] of a matrix.
+        unsafe { &*(self as *const _ as *const _) }
     }
 
     /// This transformation seen as a `Projective3`.
@@ -284,6 +285,7 @@ impl<T> Orthographic3<T> {
     /// ```
     #[inline]
     #[must_use]
+    // TODO: rename into `into_projective` to appease clippy.
     pub fn to_projective(self) -> Projective3<T> {
         Projective3::from_matrix_unchecked(self.matrix)
     }
diff --git a/src/geometry/perspective.rs b/src/geometry/perspective.rs
index 73023080..0a0e34e9 100644
--- a/src/geometry/perspective.rs
+++ b/src/geometry/perspective.rs
@@ -139,7 +139,8 @@ impl<T: RealField> Perspective3<T> {
     #[inline]
     #[must_use]
     pub fn as_projective(&self) -> &Projective3<T> {
-        unsafe { &*(self as *const _ as *const Projective3<T>) }
+        // Safety: Self and Projective3 are both #[repr(transparent)] of a matrix.
+        unsafe { &*(self as *const _ as *const _) }
     }
 
     /// This transformation seen as a `Projective3`.
diff --git a/src/geometry/point.rs b/src/geometry/point.rs
index f3c01a94..9fc8c663 100644
--- a/src/geometry/point.rs
+++ b/src/geometry/point.rs
@@ -42,7 +42,7 @@ use crate::Scalar;
 /// achieved by multiplication, e.g., `isometry * point` or `rotation * point`. Some of these transformation
 /// may have some other methods, e.g., `isometry.inverse_transform_point(&point)`. See the documentation
 /// of said transformations for details.
-#[repr(C)]
+#[repr(transparent)]
 pub struct OPoint<T, D: DimName>
 where
     DefaultAllocator: InnerAllocator<T, D>,
diff --git a/src/geometry/point_construction.rs b/src/geometry/point_construction.rs
index 581dca8d..988cc3d6 100644
--- a/src/geometry/point_construction.rs
+++ b/src/geometry/point_construction.rs
@@ -28,7 +28,7 @@ where
 {
     /// Creates a new point with uninitialized coordinates.
     #[inline]
-    pub unsafe fn new_uninitialized() -> OPoint<MaybeUninit<T>, D> {
+    pub fn new_uninitialized() -> OPoint<MaybeUninit<T>, D> {
         OPoint::from(OVector::new_uninitialized_generic(D::name(), Const::<1>))
     }
 
diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs
index 3550cbd1..bdda6e64 100755
--- a/src/geometry/quaternion.rs
+++ b/src/geometry/quaternion.rs
@@ -26,7 +26,7 @@ use crate::geometry::{Point3, Rotation};
 
 /// A quaternion. See the type alias `UnitQuaternion = Unit<Quaternion>` for a quaternion
 /// that may be used as a rotation.
-#[repr(C)]
+#[repr(transparent)]
 #[derive(Debug, Copy, Clone)]
 pub struct Quaternion<T> {
     /// This quaternion as a 4D vector of coordinates in the `[ x, y, z, w ]` storage order.
diff --git a/src/geometry/quaternion_coordinates.rs b/src/geometry/quaternion_coordinates.rs
index ba887f63..40d8ca84 100644
--- a/src/geometry/quaternion_coordinates.rs
+++ b/src/geometry/quaternion_coordinates.rs
@@ -12,13 +12,14 @@ impl<T: Scalar + SimdValue> Deref for Quaternion<T> {
 
     #[inline]
     fn deref(&self) -> &Self::Target {
-        unsafe { &*(self as *const _ as *const Self::Target) }
+        // Safety: Self and IJKW are both stored as contiguous coordinates.
+        unsafe { &*(self as *const _ as *const _) }
     }
 }
 
 impl<T: Scalar + SimdValue> DerefMut for Quaternion<T> {
     #[inline]
     fn deref_mut(&mut self) -> &mut Self::Target {
-        unsafe { &mut *(self as *mut _ as *mut Self::Target) }
+        unsafe { &mut *(self as *mut _ as *mut _) }
     }
 }
diff --git a/src/geometry/reflection.rs b/src/geometry/reflection.rs
index 79b15a30..9cd818f5 100644
--- a/src/geometry/reflection.rs
+++ b/src/geometry/reflection.rs
@@ -9,7 +9,7 @@ use simba::scalar::ComplexField;
 use crate::geometry::Point;
 
 /// A reflection wrt. a plane.
-pub struct Reflection<T, D:Dim, S> {
+pub struct Reflection<T, D: Dim, S> {
     axis: Vector<T, D, S>,
     bias: T,
 }
@@ -85,8 +85,7 @@ impl<T: ComplexField, D: Dim, S: Storage<T, D>> Reflection<T, D, S> {
         S3: StorageMut<MaybeUninit<T>, R2>,
         ShapeConstraint: DimEq<C2, D> + AreMultipliable<R2, C2, D, U1>,
     {
-        lhs.mul_to(&self.axis, work);
-        let mut work = unsafe { work.assume_init_mut() };
+        let mut work = lhs.mul_to(&self.axis, work);
 
         if !self.bias.is_zero() {
             work.add_scalar_mut(-self.bias);
@@ -107,8 +106,7 @@ impl<T: ComplexField, D: Dim, S: Storage<T, D>> Reflection<T, D, S> {
         S3: StorageMut<MaybeUninit<T>, R2>,
         ShapeConstraint: DimEq<C2, D> + AreMultipliable<R2, C2, D, U1>,
     {
-        lhs.mul_to(&self.axis, work);
-        let mut work = unsafe { work.assume_init_mut() };
+        let mut work = lhs.mul_to(&self.axis, work);
 
         if !self.bias.is_zero() {
             work.add_scalar_mut(-self.bias);
diff --git a/src/geometry/rotation.rs b/src/geometry/rotation.rs
index 04ffca71..4a74c5f2 100755
--- a/src/geometry/rotation.rs
+++ b/src/geometry/rotation.rs
@@ -54,7 +54,7 @@ use crate::geometry::Point;
 /// # Conversion
 /// * [Conversion to a matrix <span style="float:right;">`matrix`, `to_homogeneous`…</span>](#conversion-to-a-matrix)
 ///
-#[repr(C)]
+#[repr(transparent)]
 #[derive(Debug)]
 pub struct Rotation<T, const D: usize> {
     matrix: SMatrix<T, D, D>,
@@ -190,7 +190,7 @@ impl<T, const D: usize> Rotation<T, D> {
     /// A mutable reference to the underlying matrix representation of this rotation.
     #[inline]
     #[deprecated(note = "Use `.matrix_mut_unchecked()` instead.")]
-    pub unsafe fn matrix_mut(&mut self) -> &mut SMatrix<T, D, D> {
+    pub  fn matrix_mut(&mut self) -> &mut SMatrix<T, D, D> {
         &mut self.matrix
     }
 
diff --git a/src/geometry/similarity.rs b/src/geometry/similarity.rs
index 19164439..3a750656 100755
--- a/src/geometry/similarity.rs
+++ b/src/geometry/similarity.rs
@@ -22,7 +22,6 @@ use crate::base::{Const, DefaultAllocator, OMatrix, SVector, Scalar};
 use crate::geometry::{AbstractRotation, Isometry, Point, Translation};
 
 /// A similarity, i.e., an uniform scaling, followed by a rotation, followed by a translation.
-#[repr(C)]
 #[derive(Debug)]
 #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))]
 #[cfg_attr(
diff --git a/src/geometry/transform.rs b/src/geometry/transform.rs
index 14bd43ae..bf61337b 100755
--- a/src/geometry/transform.rs
+++ b/src/geometry/transform.rs
@@ -157,7 +157,7 @@ super_tcategory_impl!(
 ///
 /// It is stored as a matrix with dimensions `(D + 1, D + 1)`, e.g., it stores a 4x4 matrix for a
 /// 3D transformation.
-#[repr(C)]
+#[repr(transparent)]
 pub struct Transform<T, C: TCategory, const D: usize>
 where
     Const<D>: DimNameAdd<U1>,
diff --git a/src/geometry/translation.rs b/src/geometry/translation.rs
index 69efa4d9..ff2cf32e 100755
--- a/src/geometry/translation.rs
+++ b/src/geometry/translation.rs
@@ -21,7 +21,7 @@ use crate::base::{Const, DefaultAllocator, OMatrix, SVector, Scalar};
 use crate::geometry::Point;
 
 /// A translation.
-#[repr(C)]
+#[repr(transparent)]
 #[derive(Debug)]
 pub struct Translation<T, const D: usize> {
     /// The translation coordinates, i.e., how much is added to a point's coordinates when it is
diff --git a/src/geometry/translation_coordinates.rs b/src/geometry/translation_coordinates.rs
index 44a4c8f2..bda57f59 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 _ as *const Self::Target) }
+                unsafe { &*(self as *const _ as *const _) }
             }
         }
 
         impl<T: Scalar> DerefMut for Translation<T, $D> {
             #[inline]
             fn deref_mut(&mut self) -> &mut Self::Target {
-                unsafe { &mut *(self as *mut _ as *mut Self::Target) }
+                unsafe { &mut *(self as *mut _ as *mut _) }
             }
         }
     }
diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs
index b7cb5cd6..141034a2 100644
--- a/src/linalg/bidiagonal.rs
+++ b/src/linalg/bidiagonal.rs
@@ -130,61 +130,66 @@ where
         let mut work = Matrix::new_uninitialized_generic(nrows, Const::<1>);
 
         let upper_diagonal = nrows.value() >= ncols.value();
-        if upper_diagonal {
-            for ite in 0..dim - 1 {
+
+        // Safety: all pointers involved are valid for writes, aligned, and uninitialized.
+        unsafe {
+            if upper_diagonal {
+                for ite in 0..dim - 1 {
+                    householder::clear_column_unchecked(
+                        &mut matrix,
+                        diagonal[ite].as_mut_ptr(),
+                        ite,
+                        0,
+                        None,
+                    );
+                    householder::clear_row_unchecked(
+                        &mut matrix,
+                        off_diagonal[ite].as_mut_ptr(),
+                        &mut axis_packed,
+                        &mut work,
+                        ite,
+                        1,
+                    );
+                }
+
                 householder::clear_column_unchecked(
                     &mut matrix,
-                    diagonal[ite].as_mut_ptr(),
-                    ite,
+                    diagonal[dim - 1].as_mut_ptr(),
+                    dim - 1,
                     0,
                     None,
                 );
-                householder::clear_row_unchecked(
-                    &mut matrix,
-                    off_diagonal[ite].as_mut_ptr(),
-                    &mut axis_packed,
-                    &mut work,
-                    ite,
-                    1,
-                );
-            }
+            } else {
+                for ite in 0..dim - 1 {
+                    householder::clear_row_unchecked(
+                        &mut matrix,
+                        diagonal[ite].as_mut_ptr(),
+                        &mut axis_packed,
+                        &mut work,
+                        ite,
+                        0,
+                    );
+                    householder::clear_column_unchecked(
+                        &mut matrix,
+                        off_diagonal[ite].as_mut_ptr(),
+                        ite,
+                        1,
+                        None,
+                    );
+                }
 
-            householder::clear_column_unchecked(
-                &mut matrix,
-                diagonal[dim - 1].as_mut_ptr(),
-                dim - 1,
-                0,
-                None,
-            );
-        } else {
-            for ite in 0..dim - 1 {
                 householder::clear_row_unchecked(
                     &mut matrix,
-                    diagonal[ite].as_mut_ptr(),
+                    diagonal[dim - 1].as_mut_ptr(),
                     &mut axis_packed,
                     &mut work,
-                    ite,
+                    dim - 1,
                     0,
                 );
-                householder::clear_column_unchecked(
-                    &mut matrix,
-                    off_diagonal[ite].as_mut_ptr(),
-                    ite,
-                    1,
-                    None,
-                );
             }
-
-            householder::clear_row_unchecked(
-                &mut matrix,
-                diagonal[dim - 1].as_mut_ptr(),
-                &mut axis_packed,
-                &mut work,
-                dim - 1,
-                0,
-            );
         }
 
+        // Safety: all values have been initialized.
         unsafe {
             Bidiagonal {
                 uv: matrix,
diff --git a/src/linalg/col_piv_qr.rs b/src/linalg/col_piv_qr.rs
index 4c896587..a82f0a7b 100644
--- a/src/linalg/col_piv_qr.rs
+++ b/src/linalg/col_piv_qr.rs
@@ -86,10 +86,13 @@ where
         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: unsafe { diag.assume_init() },
+            // Safety: there's no (uninitialized) values.
+            unsafe {
+                return ColPivQR {
+                    col_piv_qr: matrix,
+                    p,
+                    diag: diag.assume_init(),
+                };
             };
         }
 
@@ -99,13 +102,19 @@ where
             matrix.swap_columns(i, col_piv);
             p.append_permutation(i, col_piv);
 
-            householder::clear_column_unchecked(&mut matrix, diag[i].as_mut_ptr(), i, 0, None);
+            // Safety: the pointer is valid for writes, aligned, and uninitialized.
+            unsafe {
+                householder::clear_column_unchecked(&mut matrix, diag[i].as_mut_ptr(), i, 0, None);
+            }
         }
 
-        ColPivQR {
-            col_piv_qr: matrix,
-            p,
-            diag: unsafe { diag.assume_init() },
+        // Safety: all values have been initialized.
+        unsafe {
+            ColPivQR {
+                col_piv_qr: matrix,
+                p,
+                diag: diag.assume_init(),
+            }
         }
     }
 
diff --git a/src/linalg/hessenberg.rs b/src/linalg/hessenberg.rs
index 6a4260bf..fc0351bf 100644
--- a/src/linalg/hessenberg.rs
+++ b/src/linalg/hessenberg.rs
@@ -111,25 +111,34 @@ where
         let mut subdiag = Matrix::new_uninitialized_generic(dim.sub(Const::<1>), Const::<1>);
 
         if dim.value() == 0 {
-            return Self {
-                hess,
-                subdiag: unsafe { subdiag.assume_init() },
-            };
+            // Safety: there's no (uninitialized) values.
+            unsafe {
+                return Self {
+                    hess,
+                    subdiag: subdiag.assume_init(),
+                };
+            }
         }
 
         for ite in 0..dim.value() - 1 {
-            householder::clear_column_unchecked(
-                &mut hess,
-                subdiag[ite].as_mut_ptr(),
-                ite,
-                1,
-                Some(work),
-            );
+            // Safety: the pointer is valid for writes, aligned, and uninitialized.
+            unsafe {
+                householder::clear_column_unchecked(
+                    &mut hess,
+                    subdiag[ite].as_mut_ptr(),
+                    ite,
+                    1,
+                    Some(work),
+                );
+            }
         }
 
-        Self {
-            hess,
-            subdiag: unsafe { subdiag.assume_init() },
+        // Safety: all values have been initialized.
+        unsafe {
+            Self {
+                hess,
+                subdiag: subdiag.assume_init(),
+            }
         }
     }
 
diff --git a/src/linalg/householder.rs b/src/linalg/householder.rs
index cb65900a..06a50d8e 100644
--- a/src/linalg/householder.rs
+++ b/src/linalg/householder.rs
@@ -45,8 +45,17 @@ pub fn reflection_axis_mut<T: ComplexField, D: Dim, S: StorageMut<T, D>>(
 
 /// Uses an householder reflection to zero out the `icol`-th column, starting with the `shift + 1`-th
 /// subdiagonal element.
+///
+/// # Safety
+/// Behavior is undefined if any of the following conditions are violated:
+///
+/// - `diag_elt` must be valid for writes.
+/// - `diag_elt` must be properly aligned.
+///
+/// Furthermore, if `diag_elt` was previously initialized, this method will leak
+/// its data.
 #[doc(hidden)]
-pub fn clear_column_unchecked<T: ComplexField, R: Dim, C: Dim>(
+pub unsafe fn clear_column_unchecked<T: ComplexField, R: Dim, C: Dim>(
     matrix: &mut OMatrix<T, R, C>,
     diag_elt: *mut T,
     icol: usize,
@@ -59,9 +68,7 @@ pub fn clear_column_unchecked<T: ComplexField, R: Dim, C: Dim>(
     let mut axis = left.rows_range_mut(icol + shift..);
 
     let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis);
-    unsafe {
-        *diag_elt = reflection_norm;
-    }
+    diag_elt.write(reflection_norm);
 
     if not_zero {
         let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
@@ -75,8 +82,17 @@ pub fn clear_column_unchecked<T: ComplexField, R: Dim, C: Dim>(
 
 /// Uses an householder reflection to zero out the `irow`-th row, ending before the `shift + 1`-th
 /// superdiagonal element.
+///
+/// # Safety
+/// Behavior is undefined if any of the following conditions are violated:
+///
+/// - `diag_elt` must be valid for writes.
+/// - `diag_elt` must be properly aligned.
+///
+/// Furthermore, if `diag_elt` was previously initialized, this method will leak
+/// its data.
 #[doc(hidden)]
-pub fn clear_row_unchecked<T: ComplexField, R: Dim, C: Dim>(
+pub unsafe fn clear_row_unchecked<T: ComplexField, R: Dim, C: Dim>(
     matrix: &mut OMatrix<T, R, C>,
     diag_elt: *mut T,
     axis_packed: &mut OVector<MaybeUninit<T>, C>,
@@ -89,13 +105,11 @@ 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_init_from(&top.columns_range(irow + shift..));
-    let mut axis = unsafe { axis.assume_init_mut() };
+    let mut axis = 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;
-    }
+    diag_elt.write(reflection_norm);
 
     if not_zero {
         let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
diff --git a/src/linalg/pow.rs b/src/linalg/pow.rs
index cb2115ad..000dc8b8 100644
--- a/src/linalg/pow.rs
+++ b/src/linalg/pow.rs
@@ -47,19 +47,24 @@ where
         // Exponentiation by squares.
         loop {
             if e % two == one {
-                self.mul_to(&multiplier, &mut buf);
+                let init_buf = self.mul_to(&multiplier, &mut buf);
+                self.copy_from(&init_buf);
+
+                // Safety: `mul_to` leaves `buf` completely initialized.
                 unsafe {
-                    self.copy_from(&buf.assume_init_ref());
+                    buf.reinitialize();
                 }
-                buf.reinitialize();
             }
 
             e /= two;
-            multiplier.mul_to(&multiplier, &mut buf);
+
+            let init_buf = multiplier.mul_to(&multiplier, &mut buf);
+            multiplier.copy_from(&init_buf);
+
+            // Safety: `mul_to` leaves `buf` completely initialized.
             unsafe {
-                multiplier.copy_from(&buf.assume_init_ref());
+                buf.reinitialize();
             }
-            buf.reinitialize();
 
             if e == zero {
                 return true;
diff --git a/src/linalg/qr.rs b/src/linalg/qr.rs
index 4b7d919c..64e14a97 100644
--- a/src/linalg/qr.rs
+++ b/src/linalg/qr.rs
@@ -94,12 +94,18 @@ where
         }
 
         for i in 0..min_nrows_ncols.value() {
-            householder::clear_column_unchecked(&mut matrix, diag[i].as_mut_ptr(), i, 0, None);
+            // Safety: the pointer is valid for writes, aligned, and uninitialized.
+            unsafe {
+                householder::clear_column_unchecked(&mut matrix, diag[i].as_mut_ptr(), i, 0, None);
+            }
         }
 
-        Self {
-            qr: matrix,
-            diag: unsafe { diag.assume_init() },
+        // Safety: all values have been initialized.
+        unsafe {
+            Self {
+                qr: matrix,
+                diag: diag.assume_init(),
+            }
         }
     }
 
diff --git a/src/proptest/mod.rs b/src/proptest/mod.rs
index a6bde56c..35410ef9 100644
--- a/src/proptest/mod.rs
+++ b/src/proptest/mod.rs
@@ -263,7 +263,7 @@ where
 }
 
 /// Same as `matrix`, but without the additional anonymous generic types
-fn matrix_<R, C, ScalarStrategy>(
+fn matrix_<R: Dim, C: Dim, ScalarStrategy>(
     value_strategy: ScalarStrategy,
     rows: DimRange<R>,
     cols: DimRange<C>,
@@ -271,8 +271,6 @@ fn matrix_<R, C, ScalarStrategy>(
 where
     ScalarStrategy: Strategy + Clone + 'static,
     ScalarStrategy::Value: Scalar,
-    R: Dim,
-    C: Dim,
     DefaultAllocator: Allocator<ScalarStrategy::Value, R, C>,
 {
     let nrows = rows.lower_bound().value()..=rows.upper_bound().value();