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 = unsafe { vec_uninit(2 * n as usize) }; + let mut eigs: Vec> = unsafe { vec_uninit2(n as usize) }; + let mut rwork: Vec> = + unsafe { vec_uninit2(2 * n as usize) }; - let mut vl: Option> = - jobvl.then(|| unsafe { vec_uninit((n * n) as usize) }); - let mut vr: Option> = - jobvr.then(|| unsafe { vec_uninit((n * n) as usize) }); + let mut vl: Option>> = + jobvl.then(|| unsafe { vec_uninit2((n * n) as usize) }); + let mut vr: Option>> = + 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 = unsafe { vec_uninit(lwork) }; + let mut work: Vec> = 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 = unsafe { vec_uninit(n as usize) }; - let mut eig_im: Vec = unsafe { vec_uninit(n as usize) }; + let mut eig_re: Vec> = unsafe { vec_uninit2(n as usize) }; + let mut eig_im: Vec> = unsafe { vec_uninit2(n as usize) }; - let mut vl: Option> = - jobvl.then(|| unsafe { vec_uninit((n * n) as usize) }); - let mut vr: Option> = - jobvr.then(|| unsafe { vec_uninit((n * n) as usize) }); + let mut vl: Option>> = + jobvl.then(|| unsafe { vec_uninit2((n * n) as usize) }); + let mut vr: Option>> = + 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 = unsafe { vec_uninit(lwork) }; + let mut work: Vec> = 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 = 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> = 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(n: usize) -> Vec { 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(n: usize) -> Vec> { + let mut v = Vec::with_capacity(n); + v.set_len(n); + v +}