Fix nalgebra-lapack.

This commit is contained in:
Sébastien Crozet 2021-08-03 17:39:45 +02:00
parent eedb860565
commit 38ac9a2f9a
9 changed files with 37 additions and 50 deletions

View File

@ -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"]

View File

@ -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]);

View File

@ -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);

View File

@ -139,10 +139,3 @@ impl ComplexHelper for Complex<f64> {
self.re
}
}
unsafe fn uninitialized_vec<T: Copy>(n: usize) -> Vec<T> {
let mut res = Vec::new();
res.reserve_exact(n);
res.set_len(n);
res
}

View File

@ -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,

View File

@ -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,

View File

@ -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<Complex<T>, 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])

View File

@ -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();

View File

@ -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,