From 38ac9a2f9abc2aa35fc3bee50050ec07ba81260b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Tue, 3 Aug 2021 17:39:45 +0200 Subject: [PATCH] Fix nalgebra-lapack. --- nalgebra-lapack/Cargo.toml | 2 +- nalgebra-lapack/src/eigen.rs | 31 +++++++++++++------------- nalgebra-lapack/src/hessenberg.rs | 6 ++--- nalgebra-lapack/src/lib.rs | 7 ------ nalgebra-lapack/src/lu.rs | 2 +- nalgebra-lapack/src/qr.rs | 6 ++--- nalgebra-lapack/src/schur.rs | 12 +++++----- nalgebra-lapack/src/svd.rs | 16 ++++++------- nalgebra-lapack/src/symmetric_eigen.rs | 5 ++--- 9 files changed, 37 insertions(+), 50 deletions(-) diff --git a/nalgebra-lapack/Cargo.toml b/nalgebra-lapack/Cargo.toml index 86825a37..0670e4b1 100644 --- a/nalgebra-lapack/Cargo.toml +++ b/nalgebra-lapack/Cargo.toml @@ -22,7 +22,7 @@ proptest-support = [ "nalgebra/proptest-support" ] arbitrary = [ "nalgebra/arbitrary" ] # For BLAS/LAPACK -default = ["netlib"] +default = ["intel-mkl"] openblas = ["lapack-src/openblas"] netlib = ["lapack-src/netlib"] accelerate = ["lapack-src/accelerate"] diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 202a1428..f6628bfe 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -77,9 +77,10 @@ where let lda = n as i32; - let mut wr = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; + // TODO: avoid the initialization? + let mut wr = Matrix::zeros_generic(nrows, Const::<1>); // TODO: Tap into the workspace. - let mut wi = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; + let mut wi = Matrix::zeros_generic(nrows, Const::<1>); let mut info = 0; let mut placeholder1 = [T::zero()]; @@ -102,14 +103,13 @@ where lapack_check!(info); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; match (left_eigenvectors, eigenvectors) { (true, true) => { - let mut vl = - unsafe { Matrix::new_uninitialized_generic(nrows, ncols).assume_init() }; - let mut vr = - unsafe { Matrix::new_uninitialized_generic(nrows, ncols).assume_init() }; + // TODO: avoid the initializations? + let mut vl = Matrix::zeros_generic(nrows, ncols); + let mut vr = Matrix::zeros_generic(nrows, ncols); T::xgeev( ljob, @@ -138,8 +138,8 @@ where } } (true, false) => { - let mut vl = - unsafe { Matrix::new_uninitialized_generic(nrows, ncols).assume_init() }; + // TODO: avoid the initialization? + let mut vl = Matrix::zeros_generic(nrows, ncols); T::xgeev( ljob, @@ -168,8 +168,8 @@ where } } (false, true) => { - let mut vr = - unsafe { Matrix::new_uninitialized_generic(nrows, ncols).assume_init() }; + // TODO: avoid the initialization? + let mut vr = Matrix::zeros_generic(nrows, ncols); T::xgeev( ljob, @@ -246,8 +246,9 @@ where let lda = n as i32; - let mut wr = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; - let mut wi = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; + // TODO: avoid the initialization? + let mut wr = Matrix::zeros_generic(nrows, Const::<1>); + let mut wi = Matrix::zeros_generic(nrows, Const::<1>); let mut info = 0; let mut placeholder1 = [T::zero()]; @@ -270,7 +271,7 @@ where lapack_panic!(info); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; T::xgeev( b'T', @@ -290,7 +291,7 @@ where ); lapack_panic!(info); - let mut res = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; + let mut res = Matrix::zeros_generic(nrows, Const::<1>); for i in 0..res.len() { res[i] = Complex::new(wr[i], wi[i]); diff --git a/nalgebra-lapack/src/hessenberg.rs b/nalgebra-lapack/src/hessenberg.rs index 0a2d125e..e05349d9 100644 --- a/nalgebra-lapack/src/hessenberg.rs +++ b/nalgebra-lapack/src/hessenberg.rs @@ -59,14 +59,12 @@ where "Unable to compute the hessenberg decomposition of an empty matrix." ); - let mut tau = unsafe { - Matrix::new_uninitialized_generic(nrows.sub(Const::<1>), Const::<1>).assume_init() - }; + let mut tau = Matrix::zeros_generic(nrows.sub(Const::<1>), Const::<1>); let mut info = 0; let lwork = T::xgehrd_work_size(n, 1, n, m.as_mut_slice(), n, tau.as_mut_slice(), &mut info); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; lapack_panic!(info); diff --git a/nalgebra-lapack/src/lib.rs b/nalgebra-lapack/src/lib.rs index 9a027772..84fa03fa 100644 --- a/nalgebra-lapack/src/lib.rs +++ b/nalgebra-lapack/src/lib.rs @@ -139,10 +139,3 @@ impl ComplexHelper for Complex { self.re } } - -unsafe fn uninitialized_vec(n: usize) -> Vec { - let mut res = Vec::new(); - res.reserve_exact(n); - res.set_len(n); - res -} diff --git a/nalgebra-lapack/src/lu.rs b/nalgebra-lapack/src/lu.rs index 5fd81771..7540c75e 100644 --- a/nalgebra-lapack/src/lu.rs +++ b/nalgebra-lapack/src/lu.rs @@ -290,7 +290,7 @@ where ); lapack_check!(info); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; T::xgetri( dim, diff --git a/nalgebra-lapack/src/qr.rs b/nalgebra-lapack/src/qr.rs index c5b5c136..895e34f3 100644 --- a/nalgebra-lapack/src/qr.rs +++ b/nalgebra-lapack/src/qr.rs @@ -56,9 +56,7 @@ where let (nrows, ncols) = m.shape_generic(); let mut info = 0; - let mut tau = unsafe { - Matrix::new_uninitialized_generic(nrows.min(ncols), Const::<1>).assume_init() - }; + let mut tau = Matrix::zeros_generic(nrows.min(ncols), Const::<1>); if nrows.value() == 0 || ncols.value() == 0 { return Self { qr: m, tau }; @@ -73,7 +71,7 @@ where &mut info, ); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; T::xgeqrf( nrows.value() as i32, diff --git a/nalgebra-lapack/src/schur.rs b/nalgebra-lapack/src/schur.rs index 82177b80..13dfc05e 100644 --- a/nalgebra-lapack/src/schur.rs +++ b/nalgebra-lapack/src/schur.rs @@ -77,9 +77,9 @@ where let mut info = 0; - let mut wr = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; - let mut wi = unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; - let mut q = unsafe { Matrix::new_uninitialized_generic(nrows, ncols).assume_init() }; + let mut wr = Matrix::zeros_generic(nrows, Const::<1>); + let mut wi = Matrix::zeros_generic(nrows, Const::<1>); + let mut q = Matrix::zeros_generic(nrows, ncols); // Placeholders: let mut bwork = [0i32]; let mut unused = 0; @@ -100,7 +100,7 @@ where ); lapack_check!(info); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; T::xgees( b'V', @@ -152,9 +152,7 @@ where where DefaultAllocator: Allocator, D>, { - let mut out = unsafe { - OVector::new_uninitialized_generic(self.t.shape_generic().0, Const::<1>).assume_init() - }; + let mut out = Matrix::zeros_generic(self.t.shape_generic().0, Const::<1>); for i in 0..out.len() { out[i] = Complex::new(self.re[i], self.im[i]) diff --git a/nalgebra-lapack/src/svd.rs b/nalgebra-lapack/src/svd.rs index aee53642..972ffa1b 100644 --- a/nalgebra-lapack/src/svd.rs +++ b/nalgebra-lapack/src/svd.rs @@ -98,9 +98,9 @@ macro_rules! svd_impl( let lda = nrows.value() as i32; - let mut u = unsafe { Matrix::new_uninitialized_generic(nrows, nrows).assume_init() }; - let mut s = unsafe { Matrix::new_uninitialized_generic(nrows.min(ncols), Const::<1>).assume_init() }; - let mut vt = unsafe { Matrix::new_uninitialized_generic(ncols, ncols).assume_init() }; + let mut u = Matrix::zeros_generic(nrows, nrows); + let mut s = Matrix::zeros_generic(nrows.min(ncols), Const::<1>); + let mut vt = Matrix::zeros_generic(ncols, ncols); let ldu = nrows.value(); let ldvt = ncols.value(); @@ -108,7 +108,7 @@ macro_rules! svd_impl( let mut work = [ 0.0 ]; let mut lwork = -1 as i32; let mut info = 0; - let mut iwork = unsafe { crate::uninitialized_vec(8 * cmp::min(nrows.value(), ncols.value())) }; + let mut iwork = vec![0; 8 * cmp::min(nrows.value(), ncols.value())]; unsafe { $lapack_func(job, nrows.value() as i32, ncols.value() as i32, m.as_mut_slice(), @@ -118,7 +118,7 @@ macro_rules! svd_impl( lapack_check!(info); lwork = work[0] as i32; - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![0.0; lwork as usize]; unsafe { $lapack_func(job, nrows.value() as i32, ncols.value() as i32, m.as_mut_slice(), @@ -253,9 +253,9 @@ macro_rules! svd_complex_impl( let min_nrows_ncols = nrows.min(ncols); - let mut u = unsafe { Matrix::new_uninitialized_generic(nrows, nrows) }; - let mut s = unsafe { Matrix::new_uninitialized_generic(min_nrows_ncols, U1) }; - let mut vt = unsafe { Matrix::new_uninitialized_generic(ncols, ncols) }; + let mut u = Matrix::zeros_generic(nrows, nrows); + let mut s = Matrix::zeros_generic(min_nrows_ncols, U1); + let mut vt = Matrix::zeros_generic(ncols, ncols); let ldu = nrows.value(); let ldvt = ncols.value(); diff --git a/nalgebra-lapack/src/symmetric_eigen.rs b/nalgebra-lapack/src/symmetric_eigen.rs index ef4ef55a..8cbe63f8 100644 --- a/nalgebra-lapack/src/symmetric_eigen.rs +++ b/nalgebra-lapack/src/symmetric_eigen.rs @@ -93,14 +93,13 @@ where let lda = n as i32; - let mut values = - unsafe { Matrix::new_uninitialized_generic(nrows, Const::<1>).assume_init() }; + let mut values = Matrix::zeros_generic(nrows, Const::<1>); let mut info = 0; let lwork = T::xsyev_work_size(jobz, b'L', n as i32, m.as_mut_slice(), lda, &mut info); lapack_check!(info); - let mut work = unsafe { crate::uninitialized_vec(lwork as usize) }; + let mut work = vec![T::zero(); lwork as usize]; T::xsyev( jobz,