From 6e5b4b0d53192f69e52d37d364be4c66f232bb6f Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Wed, 31 Aug 2022 21:36:02 +0900
Subject: [PATCH 01/10] impl AsPtr for Vec<MaybeUninit<T>>

---
 lax/src/lib.rs | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/lax/src/lib.rs b/lax/src/lib.rs
index 26b740bb..3fbe7643 100644
--- a/lax/src/lib.rs
+++ b/lax/src/lib.rs
@@ -100,6 +100,7 @@ pub use self::triangular::*;
 pub use self::tridiagonal::*;
 
 use cauchy::*;
+use std::mem::MaybeUninit;
 
 pub type Pivot = Vec<i32>;
 
@@ -150,6 +151,10 @@ impl_as_ptr!(f32, f32);
 impl_as_ptr!(f64, f64);
 impl_as_ptr!(c32, lapack_sys::__BindgenComplex<f32>);
 impl_as_ptr!(c64, lapack_sys::__BindgenComplex<f64>);
+impl_as_ptr!(MaybeUninit<f32>, f32);
+impl_as_ptr!(MaybeUninit<f64>, f64);
+impl_as_ptr!(MaybeUninit<c32>, lapack_sys::__BindgenComplex<f32>);
+impl_as_ptr!(MaybeUninit<c64>, lapack_sys::__BindgenComplex<f64>);
 
 /// Upper/Lower specification for seveal usages
 #[derive(Debug, Clone, Copy)]

From e35bdbb5a421746508465d7ac849c7ee9d12deb2 Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Wed, 31 Aug 2022 21:40:40 +0900
Subject: [PATCH 02/10] Introduce VecAssumeInit to convert Vec<MaybeUninit<T>>
 to Vec<T>

Use ManuallyDrop to implement VecAssumeInit
---
 lax/src/lib.rs | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/lax/src/lib.rs b/lax/src/lib.rs
index 3fbe7643..9b1db0a3 100644
--- a/lax/src/lib.rs
+++ b/lax/src/lib.rs
@@ -156,6 +156,30 @@ impl_as_ptr!(MaybeUninit<f64>, f64);
 impl_as_ptr!(MaybeUninit<c32>, lapack_sys::__BindgenComplex<f32>);
 impl_as_ptr!(MaybeUninit<c64>, lapack_sys::__BindgenComplex<f64>);
 
+pub(crate) trait VecAssumeInit {
+    type Target;
+    unsafe fn assume_init(self) -> Self::Target;
+}
+
+macro_rules! impl_vec_assume_init {
+    ($e:ty) => {
+        impl VecAssumeInit for Vec<MaybeUninit<$e>> {
+            type Target = Vec<$e>;
+            unsafe fn assume_init(self) -> Self::Target {
+                // FIXME use Vec::into_raw_parts instead after stablized
+                // https://doc.rust-lang.org/std/vec/struct.Vec.html#method.into_raw_parts
+                let mut me = std::mem::ManuallyDrop::new(self);
+                Vec::from_raw_parts(me.as_mut_ptr() as *mut $e, me.len(), me.capacity())
+            }
+        }
+    };
+}
+
+impl_vec_assume_init!(f32);
+impl_vec_assume_init!(f64);
+impl_vec_assume_init!(c32);
+impl_vec_assume_init!(c64);
+
 /// Upper/Lower specification for seveal usages
 #[derive(Debug, Clone, Copy)]
 #[repr(u8)]

From e6d03a542b22796196f6347346f0183822043e1c Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Wed, 31 Aug 2022 22:50:44 +0900
Subject: [PATCH 03/10] Introduce alternative `vec_uninit2`, and replace
 `vec_uninit` in eig.rs

---
 lax/src/eig.rs | 49 ++++++++++++++++++++++++++++++-------------------
 lax/src/lib.rs | 12 ++++++++++++
 2 files changed, 42 insertions(+), 19 deletions(-)

diff --git a/lax/src/eig.rs b/lax/src/eig.rs
index bf3b9f8c..ce84cd43 100644
--- a/lax/src/eig.rs
+++ b/lax/src/eig.rs
@@ -41,13 +41,14 @@ macro_rules! impl_eig_complex {
                 } else {
                     (EigenVectorFlag::Not, EigenVectorFlag::Not)
                 };
-                let mut eigs = unsafe { vec_uninit(n as usize) };
-                let mut rwork: Vec<Self::Real> = unsafe { vec_uninit(2 * n as usize) };
+                let mut eigs: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(n as usize) };
+                let mut rwork: Vec<MaybeUninit<Self::Real>> =
+                    unsafe { vec_uninit2(2 * n as usize) };
 
-                let mut vl: Option<Vec<Self>> =
-                    jobvl.then(|| unsafe { vec_uninit((n * n) as usize) });
-                let mut vr: Option<Vec<Self>> =
-                    jobvr.then(|| unsafe { vec_uninit((n * n) as usize) });
+                let mut vl: Option<Vec<MaybeUninit<Self>>> =
+                    jobvl.then(|| unsafe { vec_uninit2((n * n) as usize) });
+                let mut vr: Option<Vec<MaybeUninit<Self>>> =
+                    jobvr.then(|| unsafe { vec_uninit2((n * n) as usize) });
 
                 // calc work size
                 let mut info = 0;
@@ -74,7 +75,7 @@ macro_rules! impl_eig_complex {
 
                 // actal ev
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<Self> = unsafe { vec_uninit(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
                 let lwork = lwork as i32;
                 unsafe {
                     $ev(
@@ -96,10 +97,14 @@ macro_rules! impl_eig_complex {
                 };
                 info.as_lapack_result()?;
 
+                let eigs = unsafe { eigs.assume_init() };
+                let vr = unsafe { vr.map(|v| v.assume_init()) };
+                let mut vl = unsafe { vl.map(|v| v.assume_init()) };
+
                 // Hermite conjugate
                 if jobvl.is_calc() {
                     for c in vl.as_mut().unwrap().iter_mut() {
-                        c.im = -c.im
+                        c.im = -c.im;
                     }
                 }
 
@@ -145,13 +150,13 @@ macro_rules! impl_eig_real {
                 } else {
                     (EigenVectorFlag::Not, EigenVectorFlag::Not)
                 };
-                let mut eig_re: Vec<Self> = unsafe { vec_uninit(n as usize) };
-                let mut eig_im: Vec<Self> = unsafe { vec_uninit(n as usize) };
+                let mut eig_re: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(n as usize) };
+                let mut eig_im: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(n as usize) };
 
-                let mut vl: Option<Vec<Self>> =
-                    jobvl.then(|| unsafe { vec_uninit((n * n) as usize) });
-                let mut vr: Option<Vec<Self>> =
-                    jobvr.then(|| unsafe { vec_uninit((n * n) as usize) });
+                let mut vl: Option<Vec<MaybeUninit<Self>>> =
+                    jobvl.then(|| unsafe { vec_uninit2((n * n) as usize) });
+                let mut vr: Option<Vec<MaybeUninit<Self>>> =
+                    jobvr.then(|| unsafe { vec_uninit2((n * n) as usize) });
 
                 // calc work size
                 let mut info = 0;
@@ -178,7 +183,7 @@ macro_rules! impl_eig_real {
 
                 // actual ev
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<Self> = unsafe { vec_uninit(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
                 let lwork = lwork as i32;
                 unsafe {
                     $ev(
@@ -200,6 +205,11 @@ macro_rules! impl_eig_real {
                 };
                 info.as_lapack_result()?;
 
+                let eig_re = unsafe { eig_re.assume_init() };
+                let eig_im = unsafe { eig_im.assume_init() };
+                let vl = unsafe { vl.map(|v| v.assume_init()) };
+                let vr = unsafe { vr.map(|v| v.assume_init()) };
+
                 // reconstruct eigenvalues
                 let eigs: Vec<Self::Complex> = eig_re
                     .iter()
@@ -228,14 +238,14 @@ macro_rules! impl_eig_real {
 
                 let n = n as usize;
                 let v = vr.or(vl).unwrap();
-                let mut eigvecs = unsafe { vec_uninit(n * n) };
+                let mut eigvecs: Vec<MaybeUninit<Self::Complex>> = unsafe { vec_uninit2(n * n) };
                 let mut col = 0;
                 while col < n {
                     if eig_im[col] == 0. {
                         // The corresponding eigenvalue is real.
                         for row in 0..n {
                             let re = v[row + col * n];
-                            eigvecs[row + col * n] = Self::complex(re, 0.);
+                            eigvecs[row + col * n].write(Self::complex(re, 0.));
                         }
                         col += 1;
                     } else {
@@ -247,12 +257,13 @@ macro_rules! impl_eig_real {
                             if jobvl.is_calc() {
                                 im = -im;
                             }
-                            eigvecs[row + col * n] = Self::complex(re, im);
-                            eigvecs[row + (col + 1) * n] = Self::complex(re, -im);
+                            eigvecs[row + col * n].write(Self::complex(re, im));
+                            eigvecs[row + (col + 1) * n].write(Self::complex(re, -im));
                         }
                         col += 2;
                     }
                 }
+                let eigvecs = unsafe { eigvecs.assume_init() };
 
                 Ok((eigs, eigvecs))
             }
diff --git a/lax/src/lib.rs b/lax/src/lib.rs
index 9b1db0a3..e4fe85a0 100644
--- a/lax/src/lib.rs
+++ b/lax/src/lib.rs
@@ -281,3 +281,15 @@ unsafe fn vec_uninit<T: Sized>(n: usize) -> Vec<T> {
     v.set_len(n);
     v
 }
+
+/// Create a vector without initialization
+///
+/// Safety
+/// ------
+/// - Memory is not initialized. Do not read the memory before write.
+///
+unsafe fn vec_uninit2<T: Sized>(n: usize) -> Vec<MaybeUninit<T>> {
+    let mut v = Vec::with_capacity(n);
+    v.set_len(n);
+    v
+}

From 8c63bc7d7bcabaa70a1dc53d70026e071161d647 Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Wed, 31 Aug 2022 23:30:47 +0900
Subject: [PATCH 04/10] Use vec_uninit2 in eigh.rs

---
 lax/src/eigh.rs | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/lax/src/eigh.rs b/lax/src/eigh.rs
index a8403e90..2a112217 100644
--- a/lax/src/eigh.rs
+++ b/lax/src/eigh.rs
@@ -42,10 +42,10 @@ macro_rules! impl_eigh {
                 assert_eq!(layout.len(), layout.lda());
                 let n = layout.len();
                 let jobz = if calc_v { EigenVectorFlag::Calc } else { EigenVectorFlag::Not };
-                let mut eigs = unsafe { vec_uninit(n as usize) };
+                let mut eigs: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2(n as usize) };
 
                 $(
-                let mut $rwork_ident: Vec<Self::Real> = unsafe { vec_uninit(3 * n as usize - 2 as usize) };
+                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2(3 * n as usize - 2 as usize) };
                 )*
 
                 // calc work size
@@ -69,7 +69,7 @@ macro_rules! impl_eigh {
 
                 // actual ev
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<Self> = unsafe { vec_uninit(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
                 let lwork = lwork as i32;
                 unsafe {
                     $ev(
@@ -86,6 +86,8 @@ macro_rules! impl_eigh {
                     );
                 }
                 info.as_lapack_result()?;
+
+                let eigs = unsafe { eigs.assume_init() };
                 Ok(eigs)
             }
 
@@ -99,10 +101,10 @@ macro_rules! impl_eigh {
                 assert_eq!(layout.len(), layout.lda());
                 let n = layout.len();
                 let jobz = if calc_v { EigenVectorFlag::Calc } else { EigenVectorFlag::Not };
-                let mut eigs = unsafe { vec_uninit(n as usize) };
+                let mut eigs: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2(n as usize) };
 
                 $(
-                let mut $rwork_ident: Vec<Self::Real> = unsafe { vec_uninit(3 * n as usize - 2) };
+                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2(3 * n as usize - 2) };
                 )*
 
                 // calc work size
@@ -129,7 +131,7 @@ macro_rules! impl_eigh {
 
                 // actual evg
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<Self> = unsafe { vec_uninit(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
                 let lwork = lwork as i32;
                 unsafe {
                     $evg(
@@ -149,6 +151,7 @@ macro_rules! impl_eigh {
                     );
                 }
                 info.as_lapack_result()?;
+                let eigs = unsafe { eigs.assume_init() };
                 Ok(eigs)
             }
         }

From ee00a0dadf0db5cb65596cf1edfbbec81d5293ce Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Thu, 1 Sep 2022 11:07:26 +0900
Subject: [PATCH 05/10] Impl VecAssumeInit for any Vec<T>

---
 lax/src/lib.rs | 25 ++++++++-----------------
 1 file changed, 8 insertions(+), 17 deletions(-)

diff --git a/lax/src/lib.rs b/lax/src/lib.rs
index e4fe85a0..6ae21af2 100644
--- a/lax/src/lib.rs
+++ b/lax/src/lib.rs
@@ -161,25 +161,16 @@ pub(crate) trait VecAssumeInit {
     unsafe fn assume_init(self) -> Self::Target;
 }
 
-macro_rules! impl_vec_assume_init {
-    ($e:ty) => {
-        impl VecAssumeInit for Vec<MaybeUninit<$e>> {
-            type Target = Vec<$e>;
-            unsafe fn assume_init(self) -> Self::Target {
-                // FIXME use Vec::into_raw_parts instead after stablized
-                // https://doc.rust-lang.org/std/vec/struct.Vec.html#method.into_raw_parts
-                let mut me = std::mem::ManuallyDrop::new(self);
-                Vec::from_raw_parts(me.as_mut_ptr() as *mut $e, me.len(), me.capacity())
-            }
-        }
-    };
+impl<T> VecAssumeInit for Vec<MaybeUninit<T>> {
+    type Target = Vec<T>;
+    unsafe fn assume_init(self) -> Self::Target {
+        // FIXME use Vec::into_raw_parts instead after stablized
+        // https://doc.rust-lang.org/std/vec/struct.Vec.html#method.into_raw_parts
+        let mut me = std::mem::ManuallyDrop::new(self);
+        Vec::from_raw_parts(me.as_mut_ptr() as *mut T, me.len(), me.capacity())
+    }
 }
 
-impl_vec_assume_init!(f32);
-impl_vec_assume_init!(f64);
-impl_vec_assume_init!(c32);
-impl_vec_assume_init!(c64);
-
 /// Upper/Lower specification for seveal usages
 #[derive(Debug, Clone, Copy)]
 #[repr(u8)]

From 152e3d5ac3edd7f841e4be5b26db295fdeff8cfe Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Thu, 1 Sep 2022 11:20:31 +0900
Subject: [PATCH 06/10] Rename `transpose` to `transpose_over`, create new
 `transpose` to write uninitilized buffer

---
 lax/src/layout.rs | 75 ++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 65 insertions(+), 10 deletions(-)

diff --git a/lax/src/layout.rs b/lax/src/layout.rs
index e7ab1da4..e83923b3 100644
--- a/lax/src/layout.rs
+++ b/lax/src/layout.rs
@@ -37,7 +37,7 @@
 //! This `S` for a matrix `A` is called "leading dimension of the array A" in LAPACK document, and denoted by `lda`.
 //!
 
-use cauchy::Scalar;
+use super::*;
 
 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
 pub enum MatrixLayout {
@@ -153,7 +153,7 @@ impl MatrixLayout {
 /// ------
 /// - If size of `a` and `layout` size mismatch
 ///
-pub fn square_transpose<T: Scalar>(layout: MatrixLayout, a: &mut [T]) {
+pub fn square_transpose<T: Copy>(layout: MatrixLayout, a: &mut [T]) {
     let (m, n) = layout.size();
     let n = n as usize;
     let m = m as usize;
@@ -162,23 +162,78 @@ pub fn square_transpose<T: Scalar>(layout: MatrixLayout, a: &mut [T]) {
         for j in (i + 1)..n {
             let a_ij = a[i * n + j];
             let a_ji = a[j * m + i];
-            a[i * n + j] = a_ji.conj();
-            a[j * m + i] = a_ij.conj();
+            a[i * n + j] = a_ji;
+            a[j * m + i] = a_ij;
         }
     }
 }
 
 /// Out-place transpose for general matrix
 ///
-/// Inplace transpose of non-square matrices is hard.
-/// See also: https://en.wikipedia.org/wiki/In-place_matrix_transposition
+/// Examples
+/// ---------
+///
+/// ```rust
+/// # use lax::layout::*;
+/// let layout = MatrixLayout::C { row: 2, lda: 3 };
+/// let a = vec![1., 2., 3., 4., 5., 6.];
+/// let (l, b) = transpose(layout, &a);
+/// assert_eq!(l, MatrixLayout::F { col: 3, lda: 2 });
+/// assert_eq!(b, &[1., 4., 2., 5., 3., 6.]);
+/// ```
+///
+/// ```rust
+/// # use lax::layout::*;
+/// let layout = MatrixLayout::F { col: 2, lda: 3 };
+/// let a = vec![1., 2., 3., 4., 5., 6.];
+/// let (l, b) = transpose(layout, &a);
+/// assert_eq!(l, MatrixLayout::C { row: 3, lda: 2 });
+/// assert_eq!(b, &[1., 4., 2., 5., 3., 6.]);
+/// ```
+///
+/// Panics
+/// ------
+/// - If input array size and `layout` size mismatch
+///
+pub fn transpose<T: Copy>(layout: MatrixLayout, input: &[T]) -> (MatrixLayout, Vec<T>) {
+    let (m, n) = layout.size();
+    let transposed = layout.resized(n, m).t();
+    let m = m as usize;
+    let n = n as usize;
+    assert_eq!(input.len(), m * n);
+
+    let mut out: Vec<MaybeUninit<T>> = unsafe { vec_uninit2(m * n) };
+
+    match layout {
+        MatrixLayout::C { .. } => {
+            for i in 0..m {
+                for j in 0..n {
+                    out[j * m + i].write(input[i * n + j]);
+                }
+            }
+        }
+        MatrixLayout::F { .. } => {
+            for i in 0..m {
+                for j in 0..n {
+                    out[i * n + j].write(input[j * m + i]);
+                }
+            }
+        }
+    }
+    (transposed, unsafe { out.assume_init() })
+}
+
+/// Out-place transpose for general matrix
+///
+/// Examples
+/// ---------
 ///
 /// ```rust
 /// # use lax::layout::*;
 /// let layout = MatrixLayout::C { row: 2, lda: 3 };
 /// let a = vec![1., 2., 3., 4., 5., 6.];
 /// let mut b = vec![0.0; a.len()];
-/// let l = transpose(layout, &a, &mut b);
+/// let l = transpose_over(layout, &a, &mut b);
 /// assert_eq!(l, MatrixLayout::F { col: 3, lda: 2 });
 /// assert_eq!(b, &[1., 4., 2., 5., 3., 6.]);
 /// ```
@@ -188,16 +243,16 @@ pub fn square_transpose<T: Scalar>(layout: MatrixLayout, a: &mut [T]) {
 /// let layout = MatrixLayout::F { col: 2, lda: 3 };
 /// let a = vec![1., 2., 3., 4., 5., 6.];
 /// let mut b = vec![0.0; a.len()];
-/// let l = transpose(layout, &a, &mut b);
+/// let l = transpose_over(layout, &a, &mut b);
 /// assert_eq!(l, MatrixLayout::C { row: 3, lda: 2 });
 /// assert_eq!(b, &[1., 4., 2., 5., 3., 6.]);
 /// ```
 ///
 /// Panics
 /// ------
-/// - If size of `a` and `layout` size mismatch
+/// - If input array sizes and `layout` size mismatch
 ///
-pub fn transpose<T: Scalar>(layout: MatrixLayout, from: &[T], to: &mut [T]) -> MatrixLayout {
+pub fn transpose_over<T: Copy>(layout: MatrixLayout, from: &[T], to: &mut [T]) -> MatrixLayout {
     let (m, n) = layout.size();
     let transposed = layout.resized(n, m).t();
     let m = m as usize;

From 066939192e0cf9727abfd34d69bf47d9e83b9adc Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Thu, 1 Sep 2022 12:30:29 +0900
Subject: [PATCH 07/10] Use new `transpose` and `transpose_over`

---
 lax/src/least_squares.rs | 12 +++++++-----
 lax/src/triangular.rs    | 12 +++++++-----
 lax/src/tridiagonal.rs   |  7 ++++---
 3 files changed, 18 insertions(+), 13 deletions(-)

diff --git a/lax/src/least_squares.rs b/lax/src/least_squares.rs
index 97f9a839..527fc77d 100644
--- a/lax/src/least_squares.rs
+++ b/lax/src/least_squares.rs
@@ -68,8 +68,9 @@ macro_rules! impl_least_squares {
                 let mut a_t = None;
                 let a_layout = match a_layout {
                     MatrixLayout::C { .. } => {
-                        a_t = Some(unsafe { vec_uninit( a.len()) });
-                        transpose(a_layout, a, a_t.as_mut().unwrap())
+                        let (layout, t) = transpose(a_layout, a);
+                        a_t = Some(t);
+                        layout
                     }
                     MatrixLayout::F { .. } => a_layout,
                 };
@@ -78,8 +79,9 @@ macro_rules! impl_least_squares {
                 let mut b_t = None;
                 let b_layout = match b_layout {
                     MatrixLayout::C { .. } => {
-                        b_t = Some(unsafe { vec_uninit( b.len()) });
-                        transpose(b_layout, b, b_t.as_mut().unwrap())
+                        let (layout, t) = transpose(b_layout, b);
+                        b_t = Some(t);
+                        layout
                     }
                     MatrixLayout::F { .. } => b_layout,
                 };
@@ -149,7 +151,7 @@ macro_rules! impl_least_squares {
                 // Skip a_t -> a transpose because A has been destroyed
                 // Re-transpose b
                 if let Some(b_t) = b_t {
-                    transpose(b_layout, &b_t, b);
+                    transpose_over(b_layout, &b_t, b);
                 }
 
                 Ok(LeastSquaresOutput {
diff --git a/lax/src/triangular.rs b/lax/src/triangular.rs
index 0288d6ba..e8825758 100644
--- a/lax/src/triangular.rs
+++ b/lax/src/triangular.rs
@@ -43,8 +43,9 @@ macro_rules! impl_triangular {
                 let mut a_t = None;
                 let a_layout = match a_layout {
                     MatrixLayout::C { .. } => {
-                        a_t = Some(unsafe { vec_uninit(a.len()) });
-                        transpose(a_layout, a, a_t.as_mut().unwrap())
+                        let (layout, t) = transpose(a_layout, a);
+                        a_t = Some(t);
+                        layout
                     }
                     MatrixLayout::F { .. } => a_layout,
                 };
@@ -53,8 +54,9 @@ macro_rules! impl_triangular {
                 let mut b_t = None;
                 let b_layout = match b_layout {
                     MatrixLayout::C { .. } => {
-                        b_t = Some(unsafe { vec_uninit(b.len()) });
-                        transpose(b_layout, b, b_t.as_mut().unwrap())
+                        let (layout, t) = transpose(b_layout, b);
+                        b_t = Some(t);
+                        layout
                     }
                     MatrixLayout::F { .. } => b_layout,
                 };
@@ -82,7 +84,7 @@ macro_rules! impl_triangular {
 
                 // Re-transpose b
                 if let Some(b_t) = b_t {
-                    transpose(b_layout, &b_t, b);
+                    transpose_over(b_layout, &b_t, b);
                 }
                 Ok(())
             }
diff --git a/lax/src/tridiagonal.rs b/lax/src/tridiagonal.rs
index c80ad4b5..5d243632 100644
--- a/lax/src/tridiagonal.rs
+++ b/lax/src/tridiagonal.rs
@@ -218,8 +218,9 @@ macro_rules! impl_tridiagonal {
                 let mut b_t = None;
                 let b_layout = match b_layout {
                     MatrixLayout::C { .. } => {
-                        b_t = Some(unsafe { vec_uninit( b.len()) });
-                        transpose(b_layout, b, b_t.as_mut().unwrap())
+                        let (layout, t) = transpose(b_layout, b);
+                        b_t = Some(t);
+                        layout
                     }
                     MatrixLayout::F { .. } => b_layout,
                 };
@@ -242,7 +243,7 @@ macro_rules! impl_tridiagonal {
                 }
                 info.as_lapack_result()?;
                 if let Some(b_t) = b_t {
-                    transpose(b_layout, &b_t, b);
+                    transpose_over(b_layout, &b_t, b);
                 }
                 Ok(())
             }

From 41a3247af155a15f740b9b27c38cb240f05491df Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Thu, 1 Sep 2022 15:33:58 +0900
Subject: [PATCH 08/10] Use vec_uninit2 in qr.rs, opnorm.rs

---
 lax/src/opnorm.rs | 4 ++--
 lax/src/qr.rs     | 8 +++++---
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/lax/src/opnorm.rs b/lax/src/opnorm.rs
index ddcc2c85..ba831456 100644
--- a/lax/src/opnorm.rs
+++ b/lax/src/opnorm.rs
@@ -18,8 +18,8 @@ macro_rules! impl_opnorm {
                     MatrixLayout::F { .. } => t,
                     MatrixLayout::C { .. } => t.transpose(),
                 };
-                let mut work: Vec<Self::Real> = if matches!(t, NormType::Infinity) {
-                    unsafe { vec_uninit(m as usize) }
+                let mut work: Vec<MaybeUninit<Self::Real>> = if matches!(t, NormType::Infinity) {
+                    unsafe { vec_uninit2(m as usize) }
                 } else {
                     Vec::new()
                 };
diff --git a/lax/src/qr.rs b/lax/src/qr.rs
index 33de0372..4dc1d5a7 100644
--- a/lax/src/qr.rs
+++ b/lax/src/qr.rs
@@ -25,7 +25,7 @@ macro_rules! impl_qr {
                 let m = l.lda();
                 let n = l.len();
                 let k = m.min(n);
-                let mut tau = unsafe { vec_uninit(k as usize) };
+                let mut tau = unsafe { vec_uninit2(k as usize) };
 
                 // eval work size
                 let mut info = 0;
@@ -62,7 +62,7 @@ macro_rules! impl_qr {
 
                 // calc
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<Self> = unsafe { vec_uninit(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
                 unsafe {
                     match l {
                         MatrixLayout::F { .. } => {
@@ -93,6 +93,8 @@ macro_rules! impl_qr {
                 }
                 info.as_lapack_result()?;
 
+                let tau = unsafe { tau.assume_init() };
+
                 Ok(tau)
             }
 
@@ -134,7 +136,7 @@ macro_rules! impl_qr {
 
                 // calc
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<Self> = unsafe { vec_uninit(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
                 unsafe {
                     match l {
                         MatrixLayout::F { .. } => $gqr(

From dd425f42626375a84b17e265759c0cf1cdecc642 Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Thu, 1 Sep 2022 15:58:10 +0900
Subject: [PATCH 09/10] Use vec_uninit2 in all

---
 lax/src/least_squares.rs | 12 +++++++-----
 lax/src/lib.rs           |  2 ++
 lax/src/rcond.rs         | 11 ++++++-----
 lax/src/solve.rs         |  7 ++++---
 lax/src/solveh.rs        | 11 ++++++-----
 lax/src/svd.rs           | 15 ++++++++++-----
 lax/src/svddc.rs         | 24 ++++++++++++++----------
 lax/src/tridiagonal.rs   | 14 ++++++++------
 8 files changed, 57 insertions(+), 39 deletions(-)

diff --git a/lax/src/least_squares.rs b/lax/src/least_squares.rs
index 527fc77d..71e0dbf5 100644
--- a/lax/src/least_squares.rs
+++ b/lax/src/least_squares.rs
@@ -87,7 +87,7 @@ macro_rules! impl_least_squares {
                 };
 
                 let rcond: Self::Real = -1.;
-                let mut singular_values: Vec<Self::Real> = unsafe { vec_uninit( k as usize) };
+                let mut singular_values: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2( k as usize) };
                 let mut rank: i32 = 0;
 
                 // eval work size
@@ -120,12 +120,12 @@ macro_rules! impl_least_squares {
 
                 // calc
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<Self> = unsafe { vec_uninit(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
                 let liwork = iwork_size[0].to_usize().unwrap();
-                let mut iwork = unsafe { vec_uninit(liwork) };
+                let mut iwork: Vec<MaybeUninit<i32>> = unsafe { vec_uninit2(liwork) };
                 $(
                 let lrwork = $rwork[0].to_usize().unwrap();
-                let mut $rwork: Vec<Self::Real> = unsafe { vec_uninit(lrwork) };
+                let mut $rwork: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2(lrwork) };
                 )*
                 unsafe {
                     $gelsd(
@@ -142,12 +142,14 @@ macro_rules! impl_least_squares {
                         AsPtr::as_mut_ptr(&mut work),
                         &(lwork as i32),
                         $(AsPtr::as_mut_ptr(&mut $rwork),)*
-                        iwork.as_mut_ptr(),
+                        AsPtr::as_mut_ptr(&mut iwork),
                         &mut info,
                     );
                 }
                 info.as_lapack_result()?;
 
+                let singular_values = unsafe { singular_values.assume_init() };
+
                 // Skip a_t -> a transpose because A has been destroyed
                 // Re-transpose b
                 if let Some(b_t) = b_t {
diff --git a/lax/src/lib.rs b/lax/src/lib.rs
index 6ae21af2..bc908da6 100644
--- a/lax/src/lib.rs
+++ b/lax/src/lib.rs
@@ -147,10 +147,12 @@ macro_rules! impl_as_ptr {
         }
     };
 }
+impl_as_ptr!(i32, i32);
 impl_as_ptr!(f32, f32);
 impl_as_ptr!(f64, f64);
 impl_as_ptr!(c32, lapack_sys::__BindgenComplex<f32>);
 impl_as_ptr!(c64, lapack_sys::__BindgenComplex<f64>);
+impl_as_ptr!(MaybeUninit<i32>, i32);
 impl_as_ptr!(MaybeUninit<f32>, f32);
 impl_as_ptr!(MaybeUninit<f64>, f64);
 impl_as_ptr!(MaybeUninit<c32>, lapack_sys::__BindgenComplex<f32>);
diff --git a/lax/src/rcond.rs b/lax/src/rcond.rs
index fcd4211f..39dc4799 100644
--- a/lax/src/rcond.rs
+++ b/lax/src/rcond.rs
@@ -17,8 +17,8 @@ macro_rules! impl_rcond_real {
                 let mut rcond = Self::Real::zero();
                 let mut info = 0;
 
-                let mut work: Vec<Self> = unsafe { vec_uninit(4 * n as usize) };
-                let mut iwork = unsafe { vec_uninit(n as usize) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(4 * n as usize) };
+                let mut iwork: Vec<MaybeUninit<i32>> = unsafe { vec_uninit2(n as usize) };
                 let norm_type = match l {
                     MatrixLayout::C { .. } => NormType::Infinity,
                     MatrixLayout::F { .. } => NormType::One,
@@ -32,7 +32,7 @@ macro_rules! impl_rcond_real {
                         &anorm,
                         &mut rcond,
                         AsPtr::as_mut_ptr(&mut work),
-                        iwork.as_mut_ptr(),
+                        AsPtr::as_mut_ptr(&mut iwork),
                         &mut info,
                     )
                 };
@@ -54,8 +54,9 @@ macro_rules! impl_rcond_complex {
                 let (n, _) = l.size();
                 let mut rcond = Self::Real::zero();
                 let mut info = 0;
-                let mut work: Vec<Self> = unsafe { vec_uninit(2 * n as usize) };
-                let mut rwork: Vec<Self::Real> = unsafe { vec_uninit(2 * n as usize) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(2 * n as usize) };
+                let mut rwork: Vec<MaybeUninit<Self::Real>> =
+                    unsafe { vec_uninit2(2 * n as usize) };
                 let norm_type = match l {
                     MatrixLayout::C { .. } => NormType::Infinity,
                     MatrixLayout::F { .. } => NormType::One,
diff --git a/lax/src/solve.rs b/lax/src/solve.rs
index 9c19c874..84affe39 100644
--- a/lax/src/solve.rs
+++ b/lax/src/solve.rs
@@ -33,7 +33,7 @@ macro_rules! impl_solve {
                     return Ok(Vec::new());
                 }
                 let k = ::std::cmp::min(row, col);
-                let mut ipiv = unsafe { vec_uninit(k as usize) };
+                let mut ipiv = unsafe { vec_uninit2(k as usize) };
                 let mut info = 0;
                 unsafe {
                     $getrf(
@@ -41,11 +41,12 @@ macro_rules! impl_solve {
                         &l.len(),
                         AsPtr::as_mut_ptr(a),
                         &l.lda(),
-                        ipiv.as_mut_ptr(),
+                        AsPtr::as_mut_ptr(&mut ipiv),
                         &mut info,
                     )
                 };
                 info.as_lapack_result()?;
+                let ipiv = unsafe { ipiv.assume_init() };
                 Ok(ipiv)
             }
 
@@ -74,7 +75,7 @@ macro_rules! impl_solve {
 
                 // actual
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<Self> = unsafe { vec_uninit(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
                 unsafe {
                     $getri(
                         &l.len(),
diff --git a/lax/src/solveh.rs b/lax/src/solveh.rs
index c5259dda..09221f03 100644
--- a/lax/src/solveh.rs
+++ b/lax/src/solveh.rs
@@ -20,7 +20,7 @@ macro_rules! impl_solveh {
         impl Solveh_ for $scalar {
             fn bk(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<Pivot> {
                 let (n, _) = l.size();
-                let mut ipiv = unsafe { vec_uninit(n as usize) };
+                let mut ipiv = unsafe { vec_uninit2(n as usize) };
                 if n == 0 {
                     return Ok(Vec::new());
                 }
@@ -34,7 +34,7 @@ macro_rules! impl_solveh {
                         &n,
                         AsPtr::as_mut_ptr(a),
                         &l.lda(),
-                        ipiv.as_mut_ptr(),
+                        AsPtr::as_mut_ptr(&mut ipiv),
                         AsPtr::as_mut_ptr(&mut work_size),
                         &(-1),
                         &mut info,
@@ -44,27 +44,28 @@ macro_rules! impl_solveh {
 
                 // actual
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<Self> = unsafe { vec_uninit(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
                 unsafe {
                     $trf(
                         uplo.as_ptr(),
                         &n,
                         AsPtr::as_mut_ptr(a),
                         &l.lda(),
-                        ipiv.as_mut_ptr(),
+                        AsPtr::as_mut_ptr(&mut ipiv),
                         AsPtr::as_mut_ptr(&mut work),
                         &(lwork as i32),
                         &mut info,
                     )
                 };
                 info.as_lapack_result()?;
+                let ipiv = unsafe { ipiv.assume_init() };
                 Ok(ipiv)
             }
 
             fn invh(l: MatrixLayout, uplo: UPLO, a: &mut [Self], ipiv: &Pivot) -> Result<()> {
                 let (n, _) = l.size();
                 let mut info = 0;
-                let mut work: Vec<Self> = unsafe { vec_uninit(n as usize) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(n as usize) };
                 unsafe {
                     $tri(
                         uplo.as_ptr(),
diff --git a/lax/src/svd.rs b/lax/src/svd.rs
index 0ee56428..2968cff6 100644
--- a/lax/src/svd.rs
+++ b/lax/src/svd.rs
@@ -65,21 +65,21 @@ macro_rules! impl_svd {
 
                 let m = l.lda();
                 let mut u = match ju {
-                    FlagSVD::All => Some(unsafe { vec_uninit( (m * m) as usize) }),
+                    FlagSVD::All => Some(unsafe { vec_uninit2( (m * m) as usize) }),
                     FlagSVD::No => None,
                 };
 
                 let n = l.len();
                 let mut vt = match jvt {
-                    FlagSVD::All => Some(unsafe { vec_uninit( (n * n) as usize) }),
+                    FlagSVD::All => Some(unsafe { vec_uninit2( (n * n) as usize) }),
                     FlagSVD::No => None,
                 };
 
                 let k = std::cmp::min(m, n);
-                let mut s = unsafe { vec_uninit( k as usize) };
+                let mut s = unsafe { vec_uninit2( k as usize) };
 
                 $(
-                let mut $rwork_ident: Vec<Self::Real> = unsafe { vec_uninit( 5 * k as usize) };
+                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2( 5 * k as usize) };
                 )*
 
                 // eval work size
@@ -108,7 +108,7 @@ macro_rules! impl_svd {
 
                 // calc
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<Self> = unsafe { vec_uninit( lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2( lwork) };
                 unsafe {
                     $gesvd(
                         ju.as_ptr(),
@@ -129,6 +129,11 @@ macro_rules! impl_svd {
                     );
                 }
                 info.as_lapack_result()?;
+
+                let s = unsafe { s.assume_init() };
+                let u = u.map(|v| unsafe { v.assume_init() });
+                let vt = vt.map(|v| unsafe { v.assume_init() });
+
                 match l {
                     MatrixLayout::F { .. } => Ok(SVDOutput { s, u, vt }),
                     MatrixLayout::C { .. } => Ok(SVDOutput { s, u: vt, vt: u }),
diff --git a/lax/src/svddc.rs b/lax/src/svddc.rs
index c1198286..6a703bd6 100644
--- a/lax/src/svddc.rs
+++ b/lax/src/svddc.rs
@@ -39,7 +39,7 @@ macro_rules! impl_svddc {
                 let m = l.lda();
                 let n = l.len();
                 let k = m.min(n);
-                let mut s = unsafe { vec_uninit( k as usize) };
+                let mut s = unsafe { vec_uninit2( k as usize) };
 
                 let (u_col, vt_row) = match jobz {
                     UVTFlag::Full | UVTFlag::None => (m, n),
@@ -47,12 +47,12 @@ macro_rules! impl_svddc {
                 };
                 let (mut u, mut vt) = match jobz {
                     UVTFlag::Full => (
-                        Some(unsafe { vec_uninit( (m * m) as usize) }),
-                        Some(unsafe { vec_uninit( (n * n) as usize) }),
+                        Some(unsafe { vec_uninit2( (m * m) as usize) }),
+                        Some(unsafe { vec_uninit2( (n * n) as usize) }),
                     ),
                     UVTFlag::Some => (
-                        Some(unsafe { vec_uninit( (m * u_col) as usize) }),
-                        Some(unsafe { vec_uninit( (n * vt_row) as usize) }),
+                        Some(unsafe { vec_uninit2( (m * u_col) as usize) }),
+                        Some(unsafe { vec_uninit2( (n * vt_row) as usize) }),
                     ),
                     UVTFlag::None => (None, None),
                 };
@@ -64,12 +64,12 @@ macro_rules! impl_svddc {
                     UVTFlag::None => 7 * mn,
                     _ => std::cmp::max(5*mn*mn + 5*mn, 2*mx*mn + 2*mn*mn + mn),
                 };
-                let mut $rwork_ident: Vec<Self::Real> = unsafe { vec_uninit( lrwork) };
+                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2( lrwork) };
                 )*
 
                 // eval work size
                 let mut info = 0;
-                let mut iwork = unsafe { vec_uninit( 8 * k as usize) };
+                let mut iwork: Vec<MaybeUninit<i32>> = unsafe { vec_uninit2( 8 * k as usize) };
                 let mut work_size = [Self::zero()];
                 unsafe {
                     $gesdd(
@@ -86,7 +86,7 @@ macro_rules! impl_svddc {
                         AsPtr::as_mut_ptr(&mut work_size),
                         &(-1),
                         $(AsPtr::as_mut_ptr(&mut $rwork_ident),)*
-                        iwork.as_mut_ptr(),
+                        AsPtr::as_mut_ptr(&mut iwork),
                         &mut info,
                     );
                 }
@@ -94,7 +94,7 @@ macro_rules! impl_svddc {
 
                 // do svd
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<Self> = unsafe { vec_uninit( lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2( lwork) };
                 unsafe {
                     $gesdd(
                         jobz.as_ptr(),
@@ -110,12 +110,16 @@ macro_rules! impl_svddc {
                         AsPtr::as_mut_ptr(&mut work),
                         &(lwork as i32),
                         $(AsPtr::as_mut_ptr(&mut $rwork_ident),)*
-                        iwork.as_mut_ptr(),
+                        AsPtr::as_mut_ptr(&mut iwork),
                         &mut info,
                     );
                 }
                 info.as_lapack_result()?;
 
+                let s = unsafe { s.assume_init() };
+                let u = u.map(|v| unsafe { v.assume_init() });
+                let vt = vt.map(|v| unsafe { v.assume_init() });
+
                 match l {
                     MatrixLayout::F { .. } => Ok(SVDOutput { s, u, vt }),
                     MatrixLayout::C { .. } => Ok(SVDOutput { s, u: vt, vt: u }),
diff --git a/lax/src/tridiagonal.rs b/lax/src/tridiagonal.rs
index 5d243632..3f86bd0e 100644
--- a/lax/src/tridiagonal.rs
+++ b/lax/src/tridiagonal.rs
@@ -152,8 +152,8 @@ macro_rules! impl_tridiagonal {
         impl Tridiagonal_ for $scalar {
             fn lu_tridiagonal(mut a: Tridiagonal<Self>) -> Result<LUFactorizedTridiagonal<Self>> {
                 let (n, _) = a.l.size();
-                let mut du2 = unsafe { vec_uninit( (n - 2) as usize) };
-                let mut ipiv = unsafe { vec_uninit( n as usize) };
+                let mut du2 = unsafe { vec_uninit2( (n - 2) as usize) };
+                let mut ipiv = unsafe { vec_uninit2( n as usize) };
                 // We have to calc one-norm before LU factorization
                 let a_opnorm_one = a.opnorm_one();
                 let mut info = 0;
@@ -164,11 +164,13 @@ macro_rules! impl_tridiagonal {
                         AsPtr::as_mut_ptr(&mut a.d),
                         AsPtr::as_mut_ptr(&mut a.du),
                         AsPtr::as_mut_ptr(&mut du2),
-                        ipiv.as_mut_ptr(),
+                        AsPtr::as_mut_ptr(&mut ipiv),
                         &mut info,
                     )
                 };
                 info.as_lapack_result()?;
+                let du2 = unsafe { du2.assume_init() };
+                let ipiv = unsafe { ipiv.assume_init() };
                 Ok(LUFactorizedTridiagonal {
                     a,
                     du2,
@@ -180,9 +182,9 @@ macro_rules! impl_tridiagonal {
             fn rcond_tridiagonal(lu: &LUFactorizedTridiagonal<Self>) -> Result<Self::Real> {
                 let (n, _) = lu.a.l.size();
                 let ipiv = &lu.ipiv;
-                let mut work: Vec<Self> = unsafe { vec_uninit( 2 * n as usize) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2( 2 * n as usize) };
                 $(
-                let mut $iwork = unsafe { vec_uninit( n as usize) };
+                let mut $iwork: Vec<MaybeUninit<i32>> = unsafe { vec_uninit2( n as usize) };
                 )*
                 let mut rcond = Self::Real::zero();
                 let mut info = 0;
@@ -198,7 +200,7 @@ macro_rules! impl_tridiagonal {
                         &lu.a_opnorm_one,
                         &mut rcond,
                         AsPtr::as_mut_ptr(&mut work),
-                        $($iwork.as_mut_ptr(),)*
+                        $(AsPtr::as_mut_ptr(&mut $iwork),)*
                         &mut info,
                     );
                 }

From b5b92e3b59754745aae57991344f2e61b4ceb832 Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Thu, 1 Sep 2022 16:00:11 +0900
Subject: [PATCH 10/10] Replace vec_uninit by vec_uninit2

---
 lax/src/eig.rs           | 23 +++++++++++------------
 lax/src/eigh.rs          | 12 ++++++------
 lax/src/layout.rs        |  2 +-
 lax/src/least_squares.rs |  8 ++++----
 lax/src/lib.rs           | 14 +-------------
 lax/src/opnorm.rs        |  2 +-
 lax/src/qr.rs            |  6 +++---
 lax/src/rcond.rs         |  9 ++++-----
 lax/src/solve.rs         |  4 ++--
 lax/src/solveh.rs        |  6 +++---
 lax/src/svd.rs           | 10 +++++-----
 lax/src/svddc.rs         | 16 ++++++++--------
 lax/src/tridiagonal.rs   |  8 ++++----
 13 files changed, 53 insertions(+), 67 deletions(-)

diff --git a/lax/src/eig.rs b/lax/src/eig.rs
index ce84cd43..f11f5287 100644
--- a/lax/src/eig.rs
+++ b/lax/src/eig.rs
@@ -41,14 +41,13 @@ macro_rules! impl_eig_complex {
                 } else {
                     (EigenVectorFlag::Not, EigenVectorFlag::Not)
                 };
-                let mut eigs: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(n as usize) };
-                let mut rwork: Vec<MaybeUninit<Self::Real>> =
-                    unsafe { vec_uninit2(2 * n as usize) };
+                let mut eigs: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(n as usize) };
+                let mut rwork: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit(2 * n as usize) };
 
                 let mut vl: Option<Vec<MaybeUninit<Self>>> =
-                    jobvl.then(|| unsafe { vec_uninit2((n * n) as usize) });
+                    jobvl.then(|| unsafe { vec_uninit((n * n) as usize) });
                 let mut vr: Option<Vec<MaybeUninit<Self>>> =
-                    jobvr.then(|| unsafe { vec_uninit2((n * n) as usize) });
+                    jobvr.then(|| unsafe { vec_uninit((n * n) as usize) });
 
                 // calc work size
                 let mut info = 0;
@@ -75,7 +74,7 @@ macro_rules! impl_eig_complex {
 
                 // actal ev
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(lwork) };
                 let lwork = lwork as i32;
                 unsafe {
                     $ev(
@@ -150,13 +149,13 @@ macro_rules! impl_eig_real {
                 } else {
                     (EigenVectorFlag::Not, EigenVectorFlag::Not)
                 };
-                let mut eig_re: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(n as usize) };
-                let mut eig_im: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(n as usize) };
+                let mut eig_re: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(n as usize) };
+                let mut eig_im: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(n as usize) };
 
                 let mut vl: Option<Vec<MaybeUninit<Self>>> =
-                    jobvl.then(|| unsafe { vec_uninit2((n * n) as usize) });
+                    jobvl.then(|| unsafe { vec_uninit((n * n) as usize) });
                 let mut vr: Option<Vec<MaybeUninit<Self>>> =
-                    jobvr.then(|| unsafe { vec_uninit2((n * n) as usize) });
+                    jobvr.then(|| unsafe { vec_uninit((n * n) as usize) });
 
                 // calc work size
                 let mut info = 0;
@@ -183,7 +182,7 @@ macro_rules! impl_eig_real {
 
                 // actual ev
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(lwork) };
                 let lwork = lwork as i32;
                 unsafe {
                     $ev(
@@ -238,7 +237,7 @@ macro_rules! impl_eig_real {
 
                 let n = n as usize;
                 let v = vr.or(vl).unwrap();
-                let mut eigvecs: Vec<MaybeUninit<Self::Complex>> = unsafe { vec_uninit2(n * n) };
+                let mut eigvecs: Vec<MaybeUninit<Self::Complex>> = unsafe { vec_uninit(n * n) };
                 let mut col = 0;
                 while col < n {
                     if eig_im[col] == 0. {
diff --git a/lax/src/eigh.rs b/lax/src/eigh.rs
index 2a112217..0692f921 100644
--- a/lax/src/eigh.rs
+++ b/lax/src/eigh.rs
@@ -42,10 +42,10 @@ macro_rules! impl_eigh {
                 assert_eq!(layout.len(), layout.lda());
                 let n = layout.len();
                 let jobz = if calc_v { EigenVectorFlag::Calc } else { EigenVectorFlag::Not };
-                let mut eigs: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2(n as usize) };
+                let mut eigs: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit(n as usize) };
 
                 $(
-                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2(3 * n as usize - 2 as usize) };
+                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit(3 * n as usize - 2 as usize) };
                 )*
 
                 // calc work size
@@ -69,7 +69,7 @@ macro_rules! impl_eigh {
 
                 // actual ev
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(lwork) };
                 let lwork = lwork as i32;
                 unsafe {
                     $ev(
@@ -101,10 +101,10 @@ macro_rules! impl_eigh {
                 assert_eq!(layout.len(), layout.lda());
                 let n = layout.len();
                 let jobz = if calc_v { EigenVectorFlag::Calc } else { EigenVectorFlag::Not };
-                let mut eigs: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2(n as usize) };
+                let mut eigs: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit(n as usize) };
 
                 $(
-                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2(3 * n as usize - 2) };
+                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit(3 * n as usize - 2) };
                 )*
 
                 // calc work size
@@ -131,7 +131,7 @@ macro_rules! impl_eigh {
 
                 // actual evg
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(lwork) };
                 let lwork = lwork as i32;
                 unsafe {
                     $evg(
diff --git a/lax/src/layout.rs b/lax/src/layout.rs
index e83923b3..e695d8e7 100644
--- a/lax/src/layout.rs
+++ b/lax/src/layout.rs
@@ -202,7 +202,7 @@ pub fn transpose<T: Copy>(layout: MatrixLayout, input: &[T]) -> (MatrixLayout, V
     let n = n as usize;
     assert_eq!(input.len(), m * n);
 
-    let mut out: Vec<MaybeUninit<T>> = unsafe { vec_uninit2(m * n) };
+    let mut out: Vec<MaybeUninit<T>> = unsafe { vec_uninit(m * n) };
 
     match layout {
         MatrixLayout::C { .. } => {
diff --git a/lax/src/least_squares.rs b/lax/src/least_squares.rs
index 71e0dbf5..6be44f33 100644
--- a/lax/src/least_squares.rs
+++ b/lax/src/least_squares.rs
@@ -87,7 +87,7 @@ macro_rules! impl_least_squares {
                 };
 
                 let rcond: Self::Real = -1.;
-                let mut singular_values: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2( k as usize) };
+                let mut singular_values: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit( k as usize) };
                 let mut rank: i32 = 0;
 
                 // eval work size
@@ -120,12 +120,12 @@ macro_rules! impl_least_squares {
 
                 // calc
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(lwork) };
                 let liwork = iwork_size[0].to_usize().unwrap();
-                let mut iwork: Vec<MaybeUninit<i32>> = unsafe { vec_uninit2(liwork) };
+                let mut iwork: Vec<MaybeUninit<i32>> = unsafe { vec_uninit(liwork) };
                 $(
                 let lrwork = $rwork[0].to_usize().unwrap();
-                let mut $rwork: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2(lrwork) };
+                let mut $rwork: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit(lrwork) };
                 )*
                 unsafe {
                     $gelsd(
diff --git a/lax/src/lib.rs b/lax/src/lib.rs
index bc908da6..c8d2264d 100644
--- a/lax/src/lib.rs
+++ b/lax/src/lib.rs
@@ -269,19 +269,7 @@ impl EigenVectorFlag {
 /// ------
 /// - Memory is not initialized. Do not read the memory before write.
 ///
-unsafe fn vec_uninit<T: Sized>(n: usize) -> Vec<T> {
-    let mut v = Vec::with_capacity(n);
-    v.set_len(n);
-    v
-}
-
-/// Create a vector without initialization
-///
-/// Safety
-/// ------
-/// - Memory is not initialized. Do not read the memory before write.
-///
-unsafe fn vec_uninit2<T: Sized>(n: usize) -> Vec<MaybeUninit<T>> {
+unsafe fn vec_uninit<T: Sized>(n: usize) -> Vec<MaybeUninit<T>> {
     let mut v = Vec::with_capacity(n);
     v.set_len(n);
     v
diff --git a/lax/src/opnorm.rs b/lax/src/opnorm.rs
index ba831456..fca7704c 100644
--- a/lax/src/opnorm.rs
+++ b/lax/src/opnorm.rs
@@ -19,7 +19,7 @@ macro_rules! impl_opnorm {
                     MatrixLayout::C { .. } => t.transpose(),
                 };
                 let mut work: Vec<MaybeUninit<Self::Real>> = if matches!(t, NormType::Infinity) {
-                    unsafe { vec_uninit2(m as usize) }
+                    unsafe { vec_uninit(m as usize) }
                 } else {
                     Vec::new()
                 };
diff --git a/lax/src/qr.rs b/lax/src/qr.rs
index 4dc1d5a7..553bb606 100644
--- a/lax/src/qr.rs
+++ b/lax/src/qr.rs
@@ -25,7 +25,7 @@ macro_rules! impl_qr {
                 let m = l.lda();
                 let n = l.len();
                 let k = m.min(n);
-                let mut tau = unsafe { vec_uninit2(k as usize) };
+                let mut tau = unsafe { vec_uninit(k as usize) };
 
                 // eval work size
                 let mut info = 0;
@@ -62,7 +62,7 @@ macro_rules! impl_qr {
 
                 // calc
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(lwork) };
                 unsafe {
                     match l {
                         MatrixLayout::F { .. } => {
@@ -136,7 +136,7 @@ macro_rules! impl_qr {
 
                 // calc
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(lwork) };
                 unsafe {
                     match l {
                         MatrixLayout::F { .. } => $gqr(
diff --git a/lax/src/rcond.rs b/lax/src/rcond.rs
index 39dc4799..dfc8a941 100644
--- a/lax/src/rcond.rs
+++ b/lax/src/rcond.rs
@@ -17,8 +17,8 @@ macro_rules! impl_rcond_real {
                 let mut rcond = Self::Real::zero();
                 let mut info = 0;
 
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(4 * n as usize) };
-                let mut iwork: Vec<MaybeUninit<i32>> = unsafe { vec_uninit2(n as usize) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(4 * n as usize) };
+                let mut iwork: Vec<MaybeUninit<i32>> = unsafe { vec_uninit(n as usize) };
                 let norm_type = match l {
                     MatrixLayout::C { .. } => NormType::Infinity,
                     MatrixLayout::F { .. } => NormType::One,
@@ -54,9 +54,8 @@ macro_rules! impl_rcond_complex {
                 let (n, _) = l.size();
                 let mut rcond = Self::Real::zero();
                 let mut info = 0;
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(2 * n as usize) };
-                let mut rwork: Vec<MaybeUninit<Self::Real>> =
-                    unsafe { vec_uninit2(2 * n as usize) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(2 * n as usize) };
+                let mut rwork: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit(2 * n as usize) };
                 let norm_type = match l {
                     MatrixLayout::C { .. } => NormType::Infinity,
                     MatrixLayout::F { .. } => NormType::One,
diff --git a/lax/src/solve.rs b/lax/src/solve.rs
index 84affe39..ae76f190 100644
--- a/lax/src/solve.rs
+++ b/lax/src/solve.rs
@@ -33,7 +33,7 @@ macro_rules! impl_solve {
                     return Ok(Vec::new());
                 }
                 let k = ::std::cmp::min(row, col);
-                let mut ipiv = unsafe { vec_uninit2(k as usize) };
+                let mut ipiv = unsafe { vec_uninit(k as usize) };
                 let mut info = 0;
                 unsafe {
                     $getrf(
@@ -75,7 +75,7 @@ macro_rules! impl_solve {
 
                 // actual
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(lwork) };
                 unsafe {
                     $getri(
                         &l.len(),
diff --git a/lax/src/solveh.rs b/lax/src/solveh.rs
index 09221f03..9f65978d 100644
--- a/lax/src/solveh.rs
+++ b/lax/src/solveh.rs
@@ -20,7 +20,7 @@ macro_rules! impl_solveh {
         impl Solveh_ for $scalar {
             fn bk(l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result<Pivot> {
                 let (n, _) = l.size();
-                let mut ipiv = unsafe { vec_uninit2(n as usize) };
+                let mut ipiv = unsafe { vec_uninit(n as usize) };
                 if n == 0 {
                     return Ok(Vec::new());
                 }
@@ -44,7 +44,7 @@ macro_rules! impl_solveh {
 
                 // actual
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(lwork) };
                 unsafe {
                     $trf(
                         uplo.as_ptr(),
@@ -65,7 +65,7 @@ macro_rules! impl_solveh {
             fn invh(l: MatrixLayout, uplo: UPLO, a: &mut [Self], ipiv: &Pivot) -> Result<()> {
                 let (n, _) = l.size();
                 let mut info = 0;
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2(n as usize) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit(n as usize) };
                 unsafe {
                     $tri(
                         uplo.as_ptr(),
diff --git a/lax/src/svd.rs b/lax/src/svd.rs
index 2968cff6..8c731c7a 100644
--- a/lax/src/svd.rs
+++ b/lax/src/svd.rs
@@ -65,21 +65,21 @@ macro_rules! impl_svd {
 
                 let m = l.lda();
                 let mut u = match ju {
-                    FlagSVD::All => Some(unsafe { vec_uninit2( (m * m) as usize) }),
+                    FlagSVD::All => Some(unsafe { vec_uninit( (m * m) as usize) }),
                     FlagSVD::No => None,
                 };
 
                 let n = l.len();
                 let mut vt = match jvt {
-                    FlagSVD::All => Some(unsafe { vec_uninit2( (n * n) as usize) }),
+                    FlagSVD::All => Some(unsafe { vec_uninit( (n * n) as usize) }),
                     FlagSVD::No => None,
                 };
 
                 let k = std::cmp::min(m, n);
-                let mut s = unsafe { vec_uninit2( k as usize) };
+                let mut s = unsafe { vec_uninit( k as usize) };
 
                 $(
-                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2( 5 * k as usize) };
+                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit( 5 * k as usize) };
                 )*
 
                 // eval work size
@@ -108,7 +108,7 @@ macro_rules! impl_svd {
 
                 // calc
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2( lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit( lwork) };
                 unsafe {
                     $gesvd(
                         ju.as_ptr(),
diff --git a/lax/src/svddc.rs b/lax/src/svddc.rs
index 6a703bd6..f956d848 100644
--- a/lax/src/svddc.rs
+++ b/lax/src/svddc.rs
@@ -39,7 +39,7 @@ macro_rules! impl_svddc {
                 let m = l.lda();
                 let n = l.len();
                 let k = m.min(n);
-                let mut s = unsafe { vec_uninit2( k as usize) };
+                let mut s = unsafe { vec_uninit( k as usize) };
 
                 let (u_col, vt_row) = match jobz {
                     UVTFlag::Full | UVTFlag::None => (m, n),
@@ -47,12 +47,12 @@ macro_rules! impl_svddc {
                 };
                 let (mut u, mut vt) = match jobz {
                     UVTFlag::Full => (
-                        Some(unsafe { vec_uninit2( (m * m) as usize) }),
-                        Some(unsafe { vec_uninit2( (n * n) as usize) }),
+                        Some(unsafe { vec_uninit( (m * m) as usize) }),
+                        Some(unsafe { vec_uninit( (n * n) as usize) }),
                     ),
                     UVTFlag::Some => (
-                        Some(unsafe { vec_uninit2( (m * u_col) as usize) }),
-                        Some(unsafe { vec_uninit2( (n * vt_row) as usize) }),
+                        Some(unsafe { vec_uninit( (m * u_col) as usize) }),
+                        Some(unsafe { vec_uninit( (n * vt_row) as usize) }),
                     ),
                     UVTFlag::None => (None, None),
                 };
@@ -64,12 +64,12 @@ macro_rules! impl_svddc {
                     UVTFlag::None => 7 * mn,
                     _ => std::cmp::max(5*mn*mn + 5*mn, 2*mx*mn + 2*mn*mn + mn),
                 };
-                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit2( lrwork) };
+                let mut $rwork_ident: Vec<MaybeUninit<Self::Real>> = unsafe { vec_uninit( lrwork) };
                 )*
 
                 // eval work size
                 let mut info = 0;
-                let mut iwork: Vec<MaybeUninit<i32>> = unsafe { vec_uninit2( 8 * k as usize) };
+                let mut iwork: Vec<MaybeUninit<i32>> = unsafe { vec_uninit( 8 * k as usize) };
                 let mut work_size = [Self::zero()];
                 unsafe {
                     $gesdd(
@@ -94,7 +94,7 @@ macro_rules! impl_svddc {
 
                 // do svd
                 let lwork = work_size[0].to_usize().unwrap();
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2( lwork) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit( lwork) };
                 unsafe {
                     $gesdd(
                         jobz.as_ptr(),
diff --git a/lax/src/tridiagonal.rs b/lax/src/tridiagonal.rs
index 3f86bd0e..ef8dfdf6 100644
--- a/lax/src/tridiagonal.rs
+++ b/lax/src/tridiagonal.rs
@@ -152,8 +152,8 @@ macro_rules! impl_tridiagonal {
         impl Tridiagonal_ for $scalar {
             fn lu_tridiagonal(mut a: Tridiagonal<Self>) -> Result<LUFactorizedTridiagonal<Self>> {
                 let (n, _) = a.l.size();
-                let mut du2 = unsafe { vec_uninit2( (n - 2) as usize) };
-                let mut ipiv = unsafe { vec_uninit2( n as usize) };
+                let mut du2 = unsafe { vec_uninit( (n - 2) as usize) };
+                let mut ipiv = unsafe { vec_uninit( n as usize) };
                 // We have to calc one-norm before LU factorization
                 let a_opnorm_one = a.opnorm_one();
                 let mut info = 0;
@@ -182,9 +182,9 @@ macro_rules! impl_tridiagonal {
             fn rcond_tridiagonal(lu: &LUFactorizedTridiagonal<Self>) -> Result<Self::Real> {
                 let (n, _) = lu.a.l.size();
                 let ipiv = &lu.ipiv;
-                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit2( 2 * n as usize) };
+                let mut work: Vec<MaybeUninit<Self>> = unsafe { vec_uninit( 2 * n as usize) };
                 $(
-                let mut $iwork: Vec<MaybeUninit<i32>> = unsafe { vec_uninit2( n as usize) };
+                let mut $iwork: Vec<MaybeUninit<i32>> = unsafe { vec_uninit( n as usize) };
                 )*
                 let mut rcond = Self::Real::zero();
                 let mut info = 0;