From 36a3ac814fbb90cb28a932456b49e3beb0ba9de7 Mon Sep 17 00:00:00 2001
From: Avi Weinstock <aweinstock314@gmail.com>
Date: Thu, 26 Nov 2020 23:25:36 -0500
Subject: [PATCH] Propagate `mem::MaybeUninit` through the return types of
 `Allocator::allocate_uninitialized` and `Matrix::new_uninitialized_generic`.

Most call sites still invoke UB through `assume_init`. Said call sites instead invoke `unimplemented!()` if the `no_unsound_assume_init` feature is enabled, to make it easier to gradually fix them.

Progress towards #556.
---
 Cargo.toml                          |  1 +
 src/base/allocator.rs               |  3 +-
 src/base/blas.rs                    | 10 ++++-
 src/base/construction.rs            | 32 ++++++++++++---
 src/base/conversion.rs              | 22 ++++++----
 src/base/default_allocator.rs       | 33 ++++++++++-----
 src/base/edition.rs                 | 22 ++++++----
 src/base/matrix.rs                  | 64 +++++++++++++++++++++++------
 src/base/mod.rs                     |  2 +-
 src/base/ops.rs                     | 18 ++++----
 src/base/statistics.rs              | 10 ++++-
 src/geometry/point.rs               |  5 ++-
 src/geometry/point_construction.rs  |  5 ++-
 src/lib.rs                          |  1 +
 src/linalg/bidiagonal.rs            | 12 +++---
 src/linalg/cholesky.rs              |  8 ++--
 src/linalg/hessenberg.rs            |  6 +--
 src/linalg/permutation_sequence.rs  |  5 ++-
 src/linalg/qr.rs                    |  2 +-
 src/linalg/schur.rs                 | 12 +++---
 src/linalg/symmetric_tridiagonal.rs |  6 +--
 21 files changed, 192 insertions(+), 87 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index 314f5174..ca331014 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -36,6 +36,7 @@ compare = [ "matrixcompare-core" ]
 libm = [ "simba/libm" ]
 libm-force = [ "simba/libm_force" ]
 proptest-support = [ "proptest" ]
+no_unsound_assume_init = [ ]
 
 # This feature is only used for tests, and enables tests that require more time to run
 slow-tests = []
diff --git a/src/base/allocator.rs b/src/base/allocator.rs
index ebd55553..3632cf5d 100644
--- a/src/base/allocator.rs
+++ b/src/base/allocator.rs
@@ -1,6 +1,7 @@
 //! Abstract definition of a matrix data storage allocator.
 
 use std::any::Any;
+use std::mem;
 
 use crate::base::constraint::{SameNumberOfColumns, SameNumberOfRows, ShapeConstraint};
 use crate::base::dimension::{Dim, U1};
@@ -21,7 +22,7 @@ pub trait Allocator<N: Scalar, R: Dim, C: Dim = U1>: Any + Sized {
     type Buffer: ContiguousStorageMut<N, R, C> + Clone;
 
     /// Allocates a buffer with the given number of rows and columns without initializing its content.
-    unsafe fn allocate_uninitialized(nrows: R, ncols: C) -> Self::Buffer;
+    unsafe fn allocate_uninitialized(nrows: R, ncols: C) -> mem::MaybeUninit<Self::Buffer>;
 
     /// Allocates a buffer initialized with the content of the given iterator.
     fn allocate_from_iterator<I: IntoIterator<Item = N>>(
diff --git a/src/base/blas.rs b/src/base/blas.rs
index 761077e5..ea56c620 100644
--- a/src/base/blas.rs
+++ b/src/base/blas.rs
@@ -1328,7 +1328,10 @@ where
         ShapeConstraint: DimEq<D1, D1> + DimEq<D1, R3> + DimEq<C3, D4>,
         DefaultAllocator: Allocator<N, D1>,
     {
-        let mut work = unsafe { Vector::new_uninitialized_generic(self.data.shape().0, U1) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut work = Vector::zeros_generic(self.data.shape().0, U1);
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut work = unsafe { Vector::new_uninitialized_generic(self.data.shape().0, U1).assume_init() };
         self.quadform_tr_with_workspace(&mut work, alpha, lhs, mid, beta)
     }
 
@@ -1421,7 +1424,10 @@ where
         ShapeConstraint: DimEq<D2, R3> + DimEq<D1, C3> + AreMultipliable<C3, R3, D2, U1>,
         DefaultAllocator: Allocator<N, D2>,
     {
-        let mut work = unsafe { Vector::new_uninitialized_generic(mid.data.shape().0, U1) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut work = Vector::zeros_generic(mid.data.shape().0, U1);
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut work = unsafe { Vector::new_uninitialized_generic(mid.data.shape().0, U1).assume_init() };
         self.quadform_with_workspace(&mut work, alpha, mid, rhs, beta)
     }
 }
diff --git a/src/base/construction.rs b/src/base/construction.rs
index 8c34bf3c..722bef75 100644
--- a/src/base/construction.rs
+++ b/src/base/construction.rs
@@ -14,6 +14,7 @@ use rand::Rng;
 #[cfg(feature = "std")]
 use rand_distr::StandardNormal;
 use std::iter;
+use std::mem;
 use typenum::{self, Cmp, Greater};
 
 #[cfg(feature = "std")]
@@ -25,6 +26,16 @@ use crate::base::dimension::{Dim, DimName, Dynamic, U1, U2, U3, U4, U5, U6};
 use crate::base::storage::Storage;
 use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Scalar, Unit, Vector, VectorN};
 
+/// When "no_unsound_assume_init" is enabled, expands to `zeros_generic()` instead of `new_uninitialized_generic().assume_init()`.
+/// Intended for use in contexts where the `Scalar` type implements `num_traits::Zero`, to check whether uninitialized memory is actually performance-critical.
+#[macro_export]
+macro_rules! zero_or_uninitialized_generic {
+    ($nrows:expr, $ncols:expr) => {{
+        #[cfg(feature="no_unsound_assume_init")] { crate::base::Matrix::zeros_generic($nrows, $ncols) }
+        #[cfg(not(feature="no_unsound_assume_init"))] { crate::base::Matrix::new_uninitialized_generic($nrows, $ncols).assume_init() }
+    }}
+}
+
 /// # Generic constructors
 /// This set of matrix and vector construction functions are all generic
 /// with-regard to the matrix dimensions. They all expect to be given
@@ -38,8 +49,8 @@ where
     /// Creates a new uninitialized matrix. If the matrix has a compile-time dimension, this panics
     /// if `nrows != R::to_usize()` or `ncols != C::to_usize()`.
     #[inline]
-    pub unsafe fn new_uninitialized_generic(nrows: R, ncols: C) -> Self {
-        Self::from_data(DefaultAllocator::allocate_uninitialized(nrows, ncols))
+    pub unsafe fn new_uninitialized_generic(nrows: R, ncols: C) -> mem::MaybeUninit<Self> {
+        Self::from_uninitialized_data(DefaultAllocator::allocate_uninitialized(nrows, ncols))
     }
 
     /// Creates a matrix with all its elements set to `elem`.
@@ -88,7 +99,10 @@ where
             "Matrix init. error: the slice did not contain the right number of elements."
         );
 
-        let mut res = unsafe { Self::new_uninitialized_generic(nrows, ncols) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: Self = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { Self::new_uninitialized_generic(nrows, ncols).assume_init() };
         let mut iter = slice.iter();
 
         for i in 0..nrows.value() {
@@ -114,7 +128,10 @@ where
     where
         F: FnMut(usize, usize) -> N,
     {
-        let mut res = unsafe { Self::new_uninitialized_generic(nrows, ncols) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: Self = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { Self::new_uninitialized_generic(nrows, ncols).assume_init() };
 
         for j in 0..ncols.value() {
             for i in 0..nrows.value() {
@@ -356,7 +373,7 @@ macro_rules! impl_constructors(
     ($($Dims: ty),*; $(=> $DimIdent: ident: $DimBound: ident),*; $($gargs: expr),*; $($args: ident),*) => {
         /// Creates a new uninitialized matrix or vector.
         #[inline]
-        pub unsafe fn new_uninitialized($($args: usize),*) -> Self {
+        pub unsafe fn new_uninitialized($($args: usize),*) -> mem::MaybeUninit<Self> {
             Self::new_uninitialized_generic($($gargs),*)
         }
 
@@ -865,7 +882,10 @@ macro_rules! componentwise_constructors_impl(
             #[inline]
             pub fn new($($args: N),*) -> Self {
                 unsafe {
-                    let mut res = Self::new_uninitialized();
+                    #[cfg(feature="no_unsound_assume_init")]
+                    let mut res: Self = unimplemented!();
+                    #[cfg(not(feature="no_unsound_assume_init"))]
+                    let mut res = Self::new_uninitialized().assume_init();
                     $( *res.get_unchecked_mut(($irow, $icol)) = $args; )*
 
                     res
diff --git a/src/base/conversion.rs b/src/base/conversion.rs
index 77bf4005..97f8ff27 100644
--- a/src/base/conversion.rs
+++ b/src/base/conversion.rs
@@ -50,7 +50,10 @@ where
         let nrows2 = R2::from_usize(nrows);
         let ncols2 = C2::from_usize(ncols);
 
-        let mut res = unsafe { MatrixMN::<N2, R2, C2>::new_uninitialized_generic(nrows2, ncols2) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: MatrixMN<N2, R2, C2> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { MatrixMN::<N2, R2, C2>::new_uninitialized_generic(nrows2, ncols2).assume_init() };
         for i in 0..nrows {
             for j in 0..ncols {
                 unsafe {
@@ -73,7 +76,10 @@ where
         let nrows = R1::from_usize(nrows2);
         let ncols = C1::from_usize(ncols2);
 
-        let mut res = unsafe { Self::new_uninitialized_generic(nrows, ncols) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: Self = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { Self::new_uninitialized_generic(nrows, ncols).assume_init() };
         for i in 0..nrows2 {
             for j in 0..ncols2 {
                 unsafe {
@@ -117,9 +123,9 @@ macro_rules! impl_from_into_asref_1D(
             fn from(arr: [N; $SZ]) -> Self {
                 unsafe {
                     let mut res = Self::new_uninitialized();
-                    ptr::copy_nonoverlapping(&arr[0], res.data.ptr_mut(), $SZ);
+                    ptr::copy_nonoverlapping(&arr[0], (*res.as_mut_ptr()).data.ptr_mut(), $SZ);
 
-                    res
+                    res.assume_init()
                 }
             }
         }
@@ -184,9 +190,9 @@ macro_rules! impl_from_into_asref_2D(
             fn from(arr: [[N; $SZRows]; $SZCols]) -> Self {
                 unsafe {
                     let mut res = Self::new_uninitialized();
-                    ptr::copy_nonoverlapping(&arr[0][0], res.data.ptr_mut(), $SZRows * $SZCols);
+                    ptr::copy_nonoverlapping(&arr[0][0], (*res.as_mut_ptr()).data.ptr_mut(), $SZRows * $SZCols);
 
-                    res
+                    res.assume_init()
                 }
             }
         }
@@ -306,13 +312,13 @@ macro_rules! impl_from_into_mint_2D(
             fn from(m: mint::$MV<N>) -> Self {
                 unsafe {
                     let mut res = Self::new_uninitialized();
-                    let mut ptr = res.data.ptr_mut();
+                    let mut ptr = (*res).data.ptr_mut();
                     $(
                         ptr::copy_nonoverlapping(&m.$component.x, ptr, $SZRows);
                         ptr = ptr.offset($SZRows);
                     )*
                     let _ = ptr;
-                    res
+                    res.assume_init()
                 }
             }
         }
diff --git a/src/base/default_allocator.rs b/src/base/default_allocator.rs
index bedca471..64b60a66 100644
--- a/src/base/default_allocator.rs
+++ b/src/base/default_allocator.rs
@@ -45,9 +45,8 @@ where
     type Buffer = ArrayStorage<N, R, C>;
 
     #[inline]
-    unsafe fn allocate_uninitialized(_: R, _: C) -> Self::Buffer {
-        // TODO: Undefined behavior, see #556
-        mem::MaybeUninit::<Self::Buffer>::uninit().assume_init()
+    unsafe fn allocate_uninitialized(_: R, _: C) -> mem::MaybeUninit<Self::Buffer> {
+        mem::MaybeUninit::<Self::Buffer>::uninit()
     }
 
     #[inline]
@@ -56,7 +55,10 @@ where
         ncols: C,
         iter: I,
     ) -> Self::Buffer {
-        let mut res = unsafe { Self::allocate_uninitialized(nrows, ncols) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: Self::Buffer = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { Self::allocate_uninitialized(nrows, ncols).assume_init() };
         let mut count = 0;
 
         for (res, e) in res.iter_mut().zip(iter.into_iter()) {
@@ -80,13 +82,13 @@ impl<N: Scalar, C: Dim> Allocator<N, Dynamic, C> for DefaultAllocator {
     type Buffer = VecStorage<N, Dynamic, C>;
 
     #[inline]
-    unsafe fn allocate_uninitialized(nrows: Dynamic, ncols: C) -> Self::Buffer {
+    unsafe fn allocate_uninitialized(nrows: Dynamic, ncols: C) -> mem::MaybeUninit<Self::Buffer> {
         let mut res = Vec::new();
         let length = nrows.value() * ncols.value();
         res.reserve_exact(length);
         res.set_len(length);
 
-        VecStorage::new(nrows, ncols, res)
+        mem::MaybeUninit::new(VecStorage::new(nrows, ncols, res))
     }
 
     #[inline]
@@ -110,13 +112,13 @@ impl<N: Scalar, R: DimName> Allocator<N, R, Dynamic> for DefaultAllocator {
     type Buffer = VecStorage<N, R, Dynamic>;
 
     #[inline]
-    unsafe fn allocate_uninitialized(nrows: R, ncols: Dynamic) -> Self::Buffer {
+    unsafe fn allocate_uninitialized(nrows: R, ncols: Dynamic) -> mem::MaybeUninit<Self::Buffer> {
         let mut res = Vec::new();
         let length = nrows.value() * ncols.value();
         res.reserve_exact(length);
         res.set_len(length);
 
-        VecStorage::new(nrows, ncols, res)
+        mem::MaybeUninit::new(VecStorage::new(nrows, ncols, res))
     }
 
     #[inline]
@@ -156,7 +158,10 @@ where
         cto: CTo,
         buf: <Self as Allocator<N, RFrom, CFrom>>::Buffer,
     ) -> ArrayStorage<N, RTo, CTo> {
-        let mut res = <Self as Allocator<N, RTo, CTo>>::allocate_uninitialized(rto, cto);
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: ArrayStorage<N, RTo, CTo> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = <Self as Allocator<N, RTo, CTo>>::allocate_uninitialized(rto, cto).assume_init();
 
         let (rfrom, cfrom) = buf.shape();
 
@@ -184,7 +189,10 @@ where
         cto: CTo,
         buf: ArrayStorage<N, RFrom, CFrom>,
     ) -> VecStorage<N, Dynamic, CTo> {
-        let mut res = <Self as Allocator<N, Dynamic, CTo>>::allocate_uninitialized(rto, cto);
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: VecStorage<N, Dynamic, CTo> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = <Self as Allocator<N, Dynamic, CTo>>::allocate_uninitialized(rto, cto).assume_init();
 
         let (rfrom, cfrom) = buf.shape();
 
@@ -212,7 +220,10 @@ where
         cto: Dynamic,
         buf: ArrayStorage<N, RFrom, CFrom>,
     ) -> VecStorage<N, RTo, Dynamic> {
-        let mut res = <Self as Allocator<N, RTo, Dynamic>>::allocate_uninitialized(rto, cto);
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: VecStorage<N, RTo, Dynamic> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = <Self as Allocator<N, RTo, Dynamic>>::allocate_uninitialized(rto, cto).assume_init();
 
         let (rfrom, cfrom) = buf.shape();
 
diff --git a/src/base/edition.rs b/src/base/edition.rs
index 983bde43..40ac5e0e 100644
--- a/src/base/edition.rs
+++ b/src/base/edition.rs
@@ -54,8 +54,7 @@ impl<N: Scalar + Zero, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     {
         let irows = irows.into_iter();
         let ncols = self.data.shape().1;
-        let mut res =
-            unsafe { MatrixMN::new_uninitialized_generic(Dynamic::new(irows.len()), ncols) };
+        let mut res = unsafe { crate::zero_or_uninitialized_generic!(Dynamic::new(irows.len()), ncols) };
 
         // First, check that all the indices from irows are valid.
         // This will allow us to use unchecked access in the inner loop.
@@ -90,7 +89,7 @@ impl<N: Scalar + Zero, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         let icols = icols.into_iter();
         let nrows = self.data.shape().0;
         let mut res =
-            unsafe { MatrixMN::new_uninitialized_generic(nrows, Dynamic::new(icols.len())) };
+            unsafe { crate::zero_or_uninitialized_generic!(nrows, Dynamic::new(icols.len())) };
 
         for (destination, source) in icols.enumerate() {
             res.column_mut(destination).copy_from(&self.column(*source))
@@ -896,7 +895,10 @@ impl<N: Scalar> DMatrix<N> {
     where
         DefaultAllocator: Reallocator<N, Dynamic, Dynamic, Dynamic, Dynamic>,
     {
-        let placeholder = unsafe { Self::new_uninitialized(0, 0) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let placeholder = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let placeholder = unsafe { Self::new_uninitialized(0, 0).assume_init() };
         let old = mem::replace(self, placeholder);
         let new = old.resize(new_nrows, new_ncols, val);
         let _ = mem::replace(self, new);
@@ -919,8 +921,10 @@ where
     where
         DefaultAllocator: Reallocator<N, Dynamic, C, Dynamic, C>,
     {
-        let placeholder =
-            unsafe { Self::new_uninitialized_generic(Dynamic::new(0), self.data.shape().1) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let placeholder = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let placeholder = unsafe { Self::new_uninitialized_generic(Dynamic::new(0), self.data.shape().1).assume_init() };
         let old = mem::replace(self, placeholder);
         let new = old.resize_vertically(new_nrows, val);
         let _ = mem::replace(self, new);
@@ -943,8 +947,10 @@ where
     where
         DefaultAllocator: Reallocator<N, R, Dynamic, R, Dynamic>,
     {
-        let placeholder =
-            unsafe { Self::new_uninitialized_generic(self.data.shape().0, Dynamic::new(0)) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let placeholder = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let placeholder = unsafe { Self::new_uninitialized_generic(self.data.shape().0, Dynamic::new(0)).assume_init() };
         let old = mem::replace(self, placeholder);
         let new = old.resize_horizontally(new_ncols, val);
         let _ = mem::replace(self, new);
diff --git a/src/base/matrix.rs b/src/base/matrix.rs
index 8035d2f8..eb5aa316 100644
--- a/src/base/matrix.rs
+++ b/src/base/matrix.rs
@@ -298,6 +298,17 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         unsafe { Self::from_data_statically_unchecked(data) }
     }
 
+    /// Creates a new uninitialized matrix with the given uninitialized data
+    pub unsafe fn from_uninitialized_data(data: mem::MaybeUninit<S>) -> mem::MaybeUninit<Self> {
+        let res: Matrix<N, R, C, mem::MaybeUninit<S>> = Matrix { data, _phantoms: PhantomData };
+        let res: mem::MaybeUninit<Matrix<N, R, C, mem::MaybeUninit<S>>> = mem::MaybeUninit::new(res);
+        // safety: since we wrap the inner MaybeUninit in an outer MaybeUninit above, the fact that the `data` field is partially-uninitialized is still opaque.
+        // with s/transmute_copy/transmute/, rustc claims that `MaybeUninit<Matrix<N, R, C, MaybeUninit<S>>>` may be of a different size from `MaybeUninit<Matrix<N, R, C, S>>`
+        // but MaybeUninit's documentation says "MaybeUninit<T> is guaranteed to have the same size, alignment, and ABI as T", which implies those types should be the same size
+        let res: mem::MaybeUninit<Matrix<N, R, C, S>> = mem::transmute_copy(&res);
+        res
+    }
+
     /// The shape of this matrix returned as the tuple (number of rows, number of columns).
     ///
     /// # Examples:
@@ -496,8 +507,10 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         let nrows: SameShapeR<R, R2> = Dim::from_usize(nrows);
         let ncols: SameShapeC<C, C2> = Dim::from_usize(ncols);
 
-        let mut res: MatrixSum<N, R, C, R2, C2> =
-            unsafe { Matrix::new_uninitialized_generic(nrows, ncols) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: MatrixSum<N, R, C, R2, C2> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res: MatrixSum<N, R, C, R2, C2> = unsafe { Matrix::new_uninitialized_generic(nrows, ncols).assume_init() };
 
         // TODO: use copy_from
         for j in 0..res.ncols() {
@@ -546,7 +559,10 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         let (nrows, ncols) = self.data.shape();
 
         unsafe {
-            let mut res = Matrix::new_uninitialized_generic(ncols, nrows);
+            #[cfg(feature="no_unsound_assume_init")]
+            let mut res = unimplemented!();
+            #[cfg(not(feature="no_unsound_assume_init"))]
+            let mut res = Matrix::new_uninitialized_generic(ncols, nrows).assume_init();
             self.transpose_to(&mut res);
 
             res
@@ -564,7 +580,10 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     {
         let (nrows, ncols) = self.data.shape();
 
-        let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: MatrixMN<N2, R, C> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols).assume_init() };
 
         for j in 0..ncols.value() {
             for i in 0..nrows.value() {
@@ -608,7 +627,10 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     {
         let (nrows, ncols) = self.data.shape();
 
-        let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: MatrixMN<N2, R, C> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols).assume_init() };
 
         for j in 0..ncols.value() {
             for i in 0..nrows.value() {
@@ -635,7 +657,10 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     {
         let (nrows, ncols) = self.data.shape();
 
-        let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: MatrixMN<N3, R, C> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols).assume_init() };
 
         assert_eq!(
             (nrows.value(), ncols.value()),
@@ -676,7 +701,10 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     {
         let (nrows, ncols) = self.data.shape();
 
-        let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: MatrixMN<N4, R, C> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols).assume_init() };
 
         assert_eq!(
             (nrows.value(), ncols.value()),
@@ -1170,7 +1198,10 @@ impl<N: SimdComplexField, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S
         let (nrows, ncols) = self.data.shape();
 
         unsafe {
-            let mut res: MatrixMN<_, C, R> = Matrix::new_uninitialized_generic(ncols, nrows);
+            #[cfg(feature="no_unsound_assume_init")]
+            let mut res = unimplemented!();
+            #[cfg(not(feature="no_unsound_assume_init"))]
+            let mut res: MatrixMN<_, C, R> = Matrix::new_uninitialized_generic(ncols, nrows).assume_init();
             self.adjoint_to(&mut res);
 
             res
@@ -1311,7 +1342,10 @@ impl<N: Scalar, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
         );
 
         let dim = self.data.shape().0;
-        let mut res = unsafe { VectorN::new_uninitialized_generic(dim, U1) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: VectorN<N2, D> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { VectorN::new_uninitialized_generic(dim, U1).assume_init() };
 
         for i in 0..dim.value() {
             unsafe {
@@ -1438,7 +1472,7 @@ impl<N: Scalar + Zero, D: DimAdd<U1>, S: Storage<N, D>> Vector<N, D, S> {
     {
         let len = self.len();
         let hnrows = DimSum::<D, U1>::from_usize(len + 1);
-        let mut res = unsafe { VectorN::<N, _>::new_uninitialized_generic(hnrows, U1) };
+        let mut res: VectorN::<N, _> = unsafe { crate::zero_or_uninitialized_generic!(hnrows, U1) };
         res.generic_slice_mut((0, 0), self.data.shape())
             .copy_from(self);
         res[(len, 0)] = element;
@@ -1783,7 +1817,10 @@ impl<N: Scalar + ClosedAdd + ClosedSub + ClosedMul, R: Dim, C: Dim, S: Storage<N
                 // TODO: soooo ugly!
                 let nrows = SameShapeR::<R, R2>::from_usize(3);
                 let ncols = SameShapeC::<C, C2>::from_usize(1);
-                let mut res = Matrix::new_uninitialized_generic(nrows, ncols);
+                #[cfg(feature="no_unsound_assume_init")]
+                let mut res: MatrixCross<N, R, C, R2, C2> = unimplemented!();
+                #[cfg(not(feature="no_unsound_assume_init"))]
+                let mut res = Matrix::new_uninitialized_generic(nrows, ncols).assume_init();
 
                 let ax = self.get_unchecked((0, 0));
                 let ay = self.get_unchecked((1, 0));
@@ -1807,7 +1844,10 @@ impl<N: Scalar + ClosedAdd + ClosedSub + ClosedMul, R: Dim, C: Dim, S: Storage<N
                 // TODO: ugly!
                 let nrows = SameShapeR::<R, R2>::from_usize(1);
                 let ncols = SameShapeC::<C, C2>::from_usize(3);
-                let mut res = Matrix::new_uninitialized_generic(nrows, ncols);
+                #[cfg(feature="no_unsound_assume_init")]
+                let mut res: MatrixCross<N, R, C, R2, C2> = unimplemented!();
+                #[cfg(not(feature="no_unsound_assume_init"))]
+                let mut res = Matrix::new_uninitialized_generic(nrows, ncols).assume_init();
 
                 let ax = self.get_unchecked((0, 0));
                 let ay = self.get_unchecked((0, 1));
diff --git a/src/base/mod.rs b/src/base/mod.rs
index edea4a2d..abe2fd48 100644
--- a/src/base/mod.rs
+++ b/src/base/mod.rs
@@ -15,7 +15,7 @@ mod alias_slice;
 mod array_storage;
 mod cg;
 mod componentwise;
-mod construction;
+#[macro_use] mod construction;
 mod construction_slice;
 mod conversion;
 mod edition;
diff --git a/src/base/ops.rs b/src/base/ops.rs
index 01968b47..37d18827 100644
--- a/src/base/ops.rs
+++ b/src/base/ops.rs
@@ -331,7 +331,10 @@ macro_rules! componentwise_binop_impl(
                     let (nrows, ncols) = self.shape();
                     let nrows: SameShapeR<R1, R2> = Dim::from_usize(nrows);
                     let ncols: SameShapeC<C1, C2> = Dim::from_usize(ncols);
-                    Matrix::new_uninitialized_generic(nrows, ncols)
+                    #[cfg(feature="no_unsound_assume_init")]
+                    { unimplemented!() }
+                    #[cfg(not(feature="no_unsound_assume_init"))]
+                    { Matrix::new_uninitialized_generic(nrows, ncols).assume_init() }
                 };
 
                 self.$method_to_statically_unchecked(rhs, &mut res);
@@ -573,9 +576,7 @@ where
 
     #[inline]
     fn mul(self, rhs: &'b Matrix<N, R2, C2, SB>) -> Self::Output {
-        let mut res =
-            unsafe { Matrix::new_uninitialized_generic(self.data.shape().0, rhs.data.shape().1) };
-
+        let mut res = unsafe { crate::zero_or_uninitialized_generic!(self.data.shape().0, rhs.data.shape().1) };
         self.mul_to(rhs, &mut res);
         res
     }
@@ -684,8 +685,7 @@ where
         DefaultAllocator: Allocator<N, C1, C2>,
         ShapeConstraint: SameNumberOfRows<R1, R2>,
     {
-        let mut res =
-            unsafe { Matrix::new_uninitialized_generic(self.data.shape().1, rhs.data.shape().1) };
+        let mut res = unsafe { crate::zero_or_uninitialized_generic!(self.data.shape().1, rhs.data.shape().1) };
 
         self.tr_mul_to(rhs, &mut res);
         res
@@ -700,8 +700,7 @@ where
         DefaultAllocator: Allocator<N, C1, C2>,
         ShapeConstraint: SameNumberOfRows<R1, R2>,
     {
-        let mut res =
-            unsafe { Matrix::new_uninitialized_generic(self.data.shape().1, rhs.data.shape().1) };
+        let mut res = unsafe { crate::zero_or_uninitialized_generic!(self.data.shape().1, rhs.data.shape().1) };
 
         self.ad_mul_to(rhs, &mut res);
         res
@@ -815,8 +814,7 @@ where
         let (nrows1, ncols1) = self.data.shape();
         let (nrows2, ncols2) = rhs.data.shape();
 
-        let mut res =
-            unsafe { Matrix::new_uninitialized_generic(nrows1.mul(nrows2), ncols1.mul(ncols2)) };
+        let mut res = unsafe { crate::zero_or_uninitialized_generic!(nrows1.mul(nrows2), ncols1.mul(ncols2)) };
 
         {
             let mut data_res = res.data.ptr_mut();
diff --git a/src/base/statistics.rs b/src/base/statistics.rs
index 231f654b..6ac826bf 100644
--- a/src/base/statistics.rs
+++ b/src/base/statistics.rs
@@ -17,7 +17,10 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         DefaultAllocator: Allocator<N, U1, C>,
     {
         let ncols = self.data.shape().1;
-        let mut res = unsafe { RowVectorN::new_uninitialized_generic(U1, ncols) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: RowVectorN<N, C> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { RowVectorN::new_uninitialized_generic(U1, ncols).assume_init() };
 
         for i in 0..ncols.value() {
             // TODO: avoid bound checking of column.
@@ -42,7 +45,10 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
         DefaultAllocator: Allocator<N, C>,
     {
         let ncols = self.data.shape().1;
-        let mut res = unsafe { VectorN::new_uninitialized_generic(ncols, U1) };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: VectorN<N, C> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { VectorN::new_uninitialized_generic(ncols, U1).assume_init() };
 
         for i in 0..ncols.value() {
             // TODO: avoid bound checking of column.
diff --git a/src/geometry/point.rs b/src/geometry/point.rs
index 75410ccd..f5151d9a 100644
--- a/src/geometry/point.rs
+++ b/src/geometry/point.rs
@@ -181,7 +181,10 @@ where
         D: DimNameAdd<U1>,
         DefaultAllocator: Allocator<N, DimNameSum<D, U1>>,
     {
-        let mut res = unsafe { VectorN::<_, DimNameSum<D, U1>>::new_uninitialized() };
+        #[cfg(feature="no_unsound_assume_init")]
+        let mut res: VectorN<N, DimNameSum<D, U1>> = unimplemented!();
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        let mut res = unsafe { VectorN::<_, DimNameSum<D, U1>>::new_uninitialized().assume_init() };
         res.fixed_slice_mut::<D, U1>(0, 0).copy_from(&self.coords);
         res[(D::dim(), 0)] = N::one();
 
diff --git a/src/geometry/point_construction.rs b/src/geometry/point_construction.rs
index f567cfac..ddf453dd 100644
--- a/src/geometry/point_construction.rs
+++ b/src/geometry/point_construction.rs
@@ -24,7 +24,10 @@ where
     /// Creates a new point with uninitialized coordinates.
     #[inline]
     pub unsafe fn new_uninitialized() -> Self {
-        Self::from(VectorN::new_uninitialized())
+        #[cfg(feature="no_unsound_assume_init")]
+        { unimplemented!() }
+        #[cfg(not(feature="no_unsound_assume_init"))]
+        { Self::from(VectorN::new_uninitialized().assume_init()) }
     }
 
     /// Creates a new point with all coordinates equal to zero.
diff --git a/src/lib.rs b/src/lib.rs
index 41620a53..6fb45ef6 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -87,6 +87,7 @@ an optimized set of tools for computer graphics and physics. Those features incl
 )]
 #![cfg_attr(not(feature = "std"), no_std)]
 #![cfg_attr(all(feature = "alloc", not(feature = "std")), feature(alloc))]
+#![cfg_attr(feature = "no_unsound_assume_init", allow(unreachable_code))]
 
 #[cfg(feature = "arbitrary")]
 extern crate quickcheck;
diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs
index 3ae38432..66957251 100644
--- a/src/linalg/bidiagonal.rs
+++ b/src/linalg/bidiagonal.rs
@@ -81,11 +81,11 @@ where
             "Cannot compute the bidiagonalization of an empty matrix."
         );
 
-        let mut diagonal = unsafe { MatrixMN::new_uninitialized_generic(min_nrows_ncols, U1) };
+        let mut diagonal = unsafe { crate::zero_or_uninitialized_generic!(min_nrows_ncols, U1) };
         let mut off_diagonal =
-            unsafe { MatrixMN::new_uninitialized_generic(min_nrows_ncols.sub(U1), U1) };
-        let mut axis_packed = unsafe { MatrixMN::new_uninitialized_generic(ncols, U1) };
-        let mut work = unsafe { MatrixMN::new_uninitialized_generic(nrows, U1) };
+            unsafe { crate::zero_or_uninitialized_generic!(min_nrows_ncols.sub(U1), U1) };
+        let mut axis_packed = unsafe { crate::zero_or_uninitialized_generic!(ncols, U1) };
+        let mut work = unsafe { crate::zero_or_uninitialized_generic!(nrows, U1) };
 
         let upper_diagonal = nrows.value() >= ncols.value();
         if upper_diagonal {
@@ -239,8 +239,8 @@ where
         let min_nrows_ncols = nrows.min(ncols);
 
         let mut res = Matrix::identity_generic(min_nrows_ncols, ncols);
-        let mut work = unsafe { MatrixMN::new_uninitialized_generic(min_nrows_ncols, U1) };
-        let mut axis_packed = unsafe { MatrixMN::new_uninitialized_generic(ncols, U1) };
+        let mut work = unsafe { crate::zero_or_uninitialized_generic!(min_nrows_ncols, U1) };
+        let mut axis_packed = unsafe { crate::zero_or_uninitialized_generic!(ncols, U1) };
 
         let shift = self.axis_shift().1;
 
diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs
index bd2f9281..134c9aa3 100644
--- a/src/linalg/cholesky.rs
+++ b/src/linalg/cholesky.rs
@@ -223,9 +223,9 @@ where
 
         // loads the data into a new matrix with an additional jth row/column
         let mut chol = unsafe {
-            Matrix::new_uninitialized_generic(
+            crate::zero_or_uninitialized_generic!(
                 self.chol.data.shape().0.add(U1),
-                self.chol.data.shape().1.add(U1),
+                self.chol.data.shape().1.add(U1)
             )
         };
         chol.slice_range_mut(..j, ..j)
@@ -288,9 +288,9 @@ where
 
         // loads the data into a new matrix except for the jth row/column
         let mut chol = unsafe {
-            Matrix::new_uninitialized_generic(
+            crate::zero_or_uninitialized_generic!(
                 self.chol.data.shape().0.sub(U1),
-                self.chol.data.shape().1.sub(U1),
+                self.chol.data.shape().1.sub(U1)
             )
         };
         chol.slice_range_mut(..j, ..j)
diff --git a/src/linalg/hessenberg.rs b/src/linalg/hessenberg.rs
index beff5420..979ddfd8 100644
--- a/src/linalg/hessenberg.rs
+++ b/src/linalg/hessenberg.rs
@@ -2,7 +2,7 @@
 use serde::{Deserialize, Serialize};
 
 use crate::allocator::Allocator;
-use crate::base::{DefaultAllocator, MatrixMN, MatrixN, VectorN};
+use crate::base::{DefaultAllocator, MatrixN, VectorN};
 use crate::dimension::{DimDiff, DimSub, U1};
 use crate::storage::Storage;
 use simba::scalar::ComplexField;
@@ -48,7 +48,7 @@ where
 {
     /// Computes the Hessenberg decomposition using householder reflections.
     pub fn new(hess: MatrixN<N, D>) -> Self {
-        let mut work = unsafe { MatrixMN::new_uninitialized_generic(hess.data.shape().0, U1) };
+        let mut work = unsafe { crate::zero_or_uninitialized_generic!(hess.data.shape().0, U1) };
         Self::new_with_workspace(hess, &mut work)
     }
 
@@ -74,7 +74,7 @@ where
             "Hessenberg: invalid workspace size."
         );
 
-        let mut subdiag = unsafe { MatrixMN::new_uninitialized_generic(dim.sub(U1), U1) };
+        let mut subdiag = unsafe { crate::zero_or_uninitialized_generic!(dim.sub(U1), U1) };
 
         if dim.value() == 0 {
             return Hessenberg { hess, subdiag };
diff --git a/src/linalg/permutation_sequence.rs b/src/linalg/permutation_sequence.rs
index 47255832..75dae37d 100644
--- a/src/linalg/permutation_sequence.rs
+++ b/src/linalg/permutation_sequence.rs
@@ -70,9 +70,12 @@ where
     #[inline]
     pub fn identity_generic(dim: D) -> Self {
         unsafe {
+            #[cfg(feature="no_unsound_assume_init")]
+            unimplemented!();
+            #[cfg(not(feature="no_unsound_assume_init"))]
             Self {
                 len: 0,
-                ipiv: VectorN::new_uninitialized_generic(dim, U1),
+                ipiv: VectorN::new_uninitialized_generic(dim, U1).assume_init(),
             }
         }
     }
diff --git a/src/linalg/qr.rs b/src/linalg/qr.rs
index f404aa5a..191ccfbe 100644
--- a/src/linalg/qr.rs
+++ b/src/linalg/qr.rs
@@ -54,7 +54,7 @@ where
         let (nrows, ncols) = matrix.data.shape();
         let min_nrows_ncols = nrows.min(ncols);
 
-        let mut diag = unsafe { MatrixMN::new_uninitialized_generic(min_nrows_ncols, U1) };
+        let mut diag = unsafe { crate::zero_or_uninitialized_generic!(min_nrows_ncols, U1) };
 
         if min_nrows_ncols.value() == 0 {
             return QR { qr: matrix, diag };
diff --git a/src/linalg/schur.rs b/src/linalg/schur.rs
index 72c9b5ac..b9229047 100644
--- a/src/linalg/schur.rs
+++ b/src/linalg/schur.rs
@@ -71,7 +71,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: MatrixN<N, D>, eps: N::RealField, max_niter: usize) -> Option<Self> {
-        let mut work = unsafe { VectorN::new_uninitialized_generic(m.data.shape().0, U1) };
+        let mut work = unsafe { crate::zero_or_uninitialized_generic!(m.data.shape().0, U1) };
 
         Self::do_decompose(m, &mut work, eps, max_niter, true)
             .map(|(q, t)| Schur { q: q.unwrap(), t })
@@ -378,7 +378,7 @@ where
     ///
     /// Return `None` if some eigenvalues are complex.
     pub fn eigenvalues(&self) -> Option<VectorN<N, D>> {
-        let mut out = unsafe { VectorN::new_uninitialized_generic(self.t.data.shape().0, U1) };
+        let mut out = unsafe { crate::zero_or_uninitialized_generic!(self.t.data.shape().0, U1) };
         if Self::do_eigenvalues(&self.t, &mut out) {
             Some(out)
         } else {
@@ -392,7 +392,7 @@ where
         N: RealField,
         DefaultAllocator: Allocator<NumComplex<N>, D>,
     {
-        let mut out = unsafe { VectorN::new_uninitialized_generic(self.t.data.shape().0, U1) };
+        let mut out = unsafe { crate::zero_or_uninitialized_generic!(self.t.data.shape().0, U1) };
         Self::do_complex_eigenvalues(&self.t, &mut out);
         out
     }
@@ -503,7 +503,7 @@ where
             "Unable to compute eigenvalues of a non-square matrix."
         );
 
-        let mut work = unsafe { VectorN::new_uninitialized_generic(self.data.shape().0, U1) };
+        let mut work = unsafe { crate::zero_or_uninitialized_generic!(self.data.shape().0, U1) };
 
         // Special case for 2x2 matrices.
         if self.nrows() == 2 {
@@ -544,7 +544,7 @@ where
         DefaultAllocator: Allocator<NumComplex<N>, D>,
     {
         let dim = self.data.shape().0;
-        let mut work = unsafe { VectorN::new_uninitialized_generic(dim, U1) };
+        let mut work = unsafe { crate::zero_or_uninitialized_generic!(dim, U1) };
 
         let schur = Schur::do_decompose(
             self.clone_owned(),
@@ -554,7 +554,7 @@ where
             false,
         )
         .unwrap();
-        let mut eig = unsafe { VectorN::new_uninitialized_generic(dim, U1) };
+        let mut eig = unsafe { crate::zero_or_uninitialized_generic!(dim, U1) };
         Schur::do_complex_eigenvalues(&schur.1, &mut eig);
         eig
     }
diff --git a/src/linalg/symmetric_tridiagonal.rs b/src/linalg/symmetric_tridiagonal.rs
index e8d9fb5d..c34a60da 100644
--- a/src/linalg/symmetric_tridiagonal.rs
+++ b/src/linalg/symmetric_tridiagonal.rs
@@ -2,7 +2,7 @@
 use serde::{Deserialize, Serialize};
 
 use crate::allocator::Allocator;
-use crate::base::{DefaultAllocator, MatrixMN, MatrixN, VectorN};
+use crate::base::{DefaultAllocator, MatrixN, VectorN};
 use crate::dimension::{DimDiff, DimSub, U1};
 use crate::storage::Storage;
 use simba::scalar::ComplexField;
@@ -61,8 +61,8 @@ where
             "Unable to compute the symmetric tridiagonal decomposition of an empty matrix."
         );
 
-        let mut off_diagonal = unsafe { MatrixMN::new_uninitialized_generic(dim.sub(U1), U1) };
-        let mut p = unsafe { MatrixMN::new_uninitialized_generic(dim.sub(U1), U1) };
+        let mut off_diagonal = unsafe { crate::zero_or_uninitialized_generic!(dim.sub(U1), U1) };
+        let mut p = unsafe { crate::zero_or_uninitialized_generic!(dim.sub(U1), U1) };
 
         for i in 0..dim.value() - 1 {
             let mut m = m.rows_range_mut(i + 1..);