srenblad/nalgebra
1
0
Fork 0
forked from M-Labs/nalgebra
Code Pull Requests Activity

enabled complex eigenvalues for lapack

Browse Source
This commit is contained in:
Marc Haubenstock 2022-05-02 21:55:51 +02:00
parent b656faa233
commit 9f8cfa207a
1 changed files with 129 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

157
nalgebra-lapack/src/eigen.rs
View File

@ -69,8 +69,8 @@ where
"Unable to compute the eigenvalue decomposition of a non-square matrix." "Unable to compute the eigenvalue decomposition of a non-square matrix."
); );
let ljob = if left_eigenvectors { b'V' } else { b'T' }; let ljob = if left_eigenvectors { b'V' } else { b'N' };
let rjob = if eigenvectors { b'V' } else { b'T' }; let rjob = if eigenvectors { b'V' } else { b'N' };
let (nrows, ncols) = m.shape_generic(); let (nrows, ncols) = m.shape_generic();
let n = nrows.value(); let n = nrows.value();
@ -232,22 +232,27 @@ where
/// The complex eigenvalues of the given matrix. /// The complex eigenvalues of the given matrix.
/// ///
/// Panics if the eigenvalue computation does not converge. /// Panics if the eigenvalue computation does not converge.
pub fn complex_eigenvalues(mut m: OMatrix<T, D, D>) -> OVector<Complex<T>, D> pub fn complex_eigenvalues(mut m: OMatrix<T, D, D>, left_eigenvectors: bool, eigenvectors: bool)
-> (OVector<Complex<T>, D>, Option<OMatrix<T, D, D>>, Option<OMatrix<T, D, D>>)
where where
DefaultAllocator: Allocator<Complex<T>, D>, DefaultAllocator: Allocator<Complex<T>, D> + Allocator<Complex<T>, D, D>,
{ {
assert!( assert!(
m.is_square(), m.is_square(),
"Unable to compute the eigenvalue decomposition of a non-square matrix." "Unable to compute the eigenvalue decomposition of a non-square matrix."
); );
let nrows = m.shape_generic().0; let ljob = if left_eigenvectors { b'V' } else { b'N' };
let rjob = if eigenvectors { b'V' } else { b'N' };
let (nrows, ncols) = m.shape_generic();
let n = nrows.value(); let n = nrows.value();
let lda = n as i32; let lda = n as i32;
// TODO: avoid the initialization? // TODO: avoid the initialization?
let mut wr = Matrix::zeros_generic(nrows, Const::<1>); let mut wr = Matrix::zeros_generic(nrows, Const::<1>);
// TODO: Tap into the workspace.
let mut wi = Matrix::zeros_generic(nrows, Const::<1>); let mut wi = Matrix::zeros_generic(nrows, Const::<1>);
let mut info = 0; let mut info = 0;
@ -255,8 +260,8 @@ where
let mut placeholder2 = [T::zero()]; let mut placeholder2 = [T::zero()];
let lwork = T::xgeev_work_size( let lwork = T::xgeev_work_size(
b'T', ljob,
b'T', rjob,
n as i32, n as i32,
m.as_mut_slice(), m.as_mut_slice(),
lda, lda,
@ -273,31 +278,127 @@ where
let mut work = vec![T::zero(); lwork as usize]; let mut work = vec![T::zero(); lwork as usize];
T::xgeev( match (left_eigenvectors, eigenvectors) {
b'T', (true, true) => {
b'T', // TODO: avoid the initializations?
n as i32, let mut vl = Matrix::zeros_generic(nrows, ncols);
m.as_mut_slice(), let mut vr = Matrix::zeros_generic(nrows, ncols);
lda,
wr.as_mut_slice(),
wi.as_mut_slice(),
&mut placeholder1,
1 as i32,
&mut placeholder2,
1 as i32,
&mut work,
lwork,
&mut info,
);
lapack_panic!(info);
let mut res = Matrix::zeros_generic(nrows, Const::<1>); T::xgeev(
ljob,
rjob,
n as i32,
m.as_mut_slice(),
lda,
wr.as_mut_slice(),
wi.as_mut_slice(),
&mut vl.as_mut_slice(),
n as i32,
&mut vr.as_mut_slice(),
n as i32,
&mut work,
lwork,
&mut info,
);
lapack_panic!(info);
for i in 0..res.len() { let mut res = Matrix::zeros_generic(nrows, Const::<1>);
res[i] = Complex::new(wr[i].clone(), wi[i].clone());
for i in 0..res.len() {
res[i] = Complex::new(wr[i].clone(), wi[i].clone());
}
return (res, Some(vl), Some(vr))
}
(true, false) => {
// TODO: avoid the initialization?
let mut vl = Matrix::zeros_generic(nrows, ncols);
T::xgeev(
ljob,
rjob,
n as i32,
m.as_mut_slice(),
lda,
wr.as_mut_slice(),
wi.as_mut_slice(),
&mut vl.as_mut_slice(),
n as i32,
&mut placeholder2,
1 as i32,
&mut work,
lwork,
&mut info,
);
lapack_panic!(info);
let mut res = Matrix::zeros_generic(nrows, Const::<1>);
for i in 0..res.len() {
res[i] = Complex::new(wr[i].clone(), wi[i].clone());
}
return (res, Some(vl), None)
}
(false, true) => {
// TODO: avoid the initialization?
let mut vr = Matrix::zeros_generic(nrows, ncols);
T::xgeev(
ljob,
rjob,
n as i32,
m.as_mut_slice(),
lda,
wr.as_mut_slice(),
wi.as_mut_slice(),
&mut placeholder1,
1 as i32,
&mut vr.as_mut_slice(),
n as i32,
&mut work,
lwork,
&mut info,
);
lapack_panic!(info);
let mut res = Matrix::zeros_generic(nrows, Const::<1>);
for i in 0..res.len() {
res[i] = Complex::new(wr[i].clone(), wi[i].clone());
}
return (res, None, Some(vr))
}
(false, false) => {
T::xgeev(
ljob,
rjob,
n as i32,
m.as_mut_slice(),
lda,
wr.as_mut_slice(),
wi.as_mut_slice(),
&mut placeholder1,
1 as i32,
&mut placeholder2,
1 as i32,
&mut work,
lwork,
&mut info,
);
lapack_panic!(info);
let mut res = Matrix::zeros_generic(nrows, Const::<1>);
for i in 0..res.len() {
res[i] = Complex::new(wr[i].clone(), wi[i].clone());
}
return (res, None, None)
}
} }
res
} }
/// The determinant of the decomposed matrix. /// The determinant of the decomposed matrix.