nalgebra-lapack: merge both Eigen decomposition function into a single one.
This commit is contained in:
Sébastien Crozet
parent
a0412c39f2
commit
3d52327f82
|
|
@ -13,7 +13,7 @@ use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar};
|
||||||
|
|
||||||
use lapack;
|
use lapack;
|
||||||
|
|
||||||
/// Eigendecomposition of a real square matrix with real eigenvalues.
|
/// Eigendecomposition of a real square matrix with real or complex eigenvalues.
|
||||||
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
|
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
|
||||||
#[cfg_attr(
|
#[cfg_attr(
|
||||||
feature = "serde-serialize",
|
feature = "serde-serialize",
|
||||||
|
|
@ -36,8 +36,10 @@ pub struct Eigen<T: Scalar, D: Dim>
|
||||||
where
|
where
|
||||||
DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
|
DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
|
||||||
{
|
{
|
||||||
/// The eigenvalues of the decomposed matrix.
|
/// The real parts of eigenvalues of the decomposed matrix.
|
||||||
pub eigenvalues: OVector<T, D>,
|
pub eigenvalues_re: OVector<T, D>,
|
||||||
|
/// The imaginary parts of the eigenvalues of the decomposed matrix.
|
||||||
|
pub eigenvalues_im: OVector<T, D>,
|
||||||
/// The (right) eigenvectors of the decomposed matrix.
|
/// The (right) eigenvectors of the decomposed matrix.
|
||||||
pub eigenvectors: Option<OMatrix<T, D, D>>,
|
pub eigenvectors: Option<OMatrix<T, D, D>>,
|
||||||
/// The left eigenvectors of the decomposed matrix.
|
/// The left eigenvectors of the decomposed matrix.
|
||||||
|
|
@ -104,169 +106,27 @@ where
|
||||||
lapack_check!(info);
|
lapack_check!(info);
|
||||||
|
|
||||||
let mut work = vec![T::zero(); lwork as usize];
|
let mut work = vec![T::zero(); lwork as usize];
|
||||||
|
let mut vl = if left_eigenvectors {
|
||||||
|
Some(Matrix::zeros_generic(nrows, ncols))
|
||||||
|
} else {
|
||||||
|
None
|
||||||
|
};
|
||||||
|
let mut vr = if eigenvectors {
|
||||||
|
Some(Matrix::zeros_generic(nrows, ncols))
|
||||||
|
} else {
|
||||||
|
None
|
||||||
|
};
|
||||||
|
|
||||||
match (left_eigenvectors, eigenvectors) {
|
let vl_ref = vl
|
||||||
(true, true) => {
|
.as_mut()
|
||||||
// TODO: avoid the initializations?
|
.map(|m| m.as_mut_slice())
|
||||||
let mut vl = Matrix::zeros_generic(nrows, ncols);
|
.unwrap_or(&mut placeholder1);
|
||||||
let mut vr = Matrix::zeros_generic(nrows, ncols);
|
let vr_ref = vr
|
||||||
|
.as_mut()
|
||||||
|
.map(|m| m.as_mut_slice())
|
||||||
|
.unwrap_or(&mut placeholder2);
|
||||||
|
|
||||||
T::xgeev(
|
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_check!(info);
|
|
||||||
|
|
||||||
if wi.iter().all(|e| e.is_zero()) {
|
|
||||||
return Some(Self {
|
|
||||||
eigenvalues: wr,
|
|
||||||
left_eigenvectors: Some(vl),
|
|
||||||
eigenvectors: 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_check!(info);
|
|
||||||
|
|
||||||
if wi.iter().all(|e| e.is_zero()) {
|
|
||||||
return Some(Self {
|
|
||||||
eigenvalues: wr,
|
|
||||||
left_eigenvectors: Some(vl),
|
|
||||||
eigenvectors: 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_check!(info);
|
|
||||||
|
|
||||||
if wi.iter().all(|e| e.is_zero()) {
|
|
||||||
return Some(Self {
|
|
||||||
eigenvalues: wr,
|
|
||||||
left_eigenvectors: None,
|
|
||||||
eigenvectors: 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_check!(info);
|
|
||||||
|
|
||||||
if wi.iter().all(|e| e.is_zero()) {
|
|
||||||
return Some(Self {
|
|
||||||
eigenvalues: wr,
|
|
||||||
left_eigenvectors: None,
|
|
||||||
eigenvectors: None,
|
|
||||||
});
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
None
|
|
||||||
}
|
|
||||||
|
|
||||||
/// The complex eigen decomposition of the given matrix.
|
|
||||||
///
|
|
||||||
/// Panics if the eigenvalue computation does not converge.
|
|
||||||
pub fn complex_eigen_decomposition(
|
|
||||||
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
|
|
||||||
DefaultAllocator: Allocator<Complex<T>, D> + Allocator<Complex<T>, D, D>,
|
|
||||||
{
|
|
||||||
assert!(
|
|
||||||
m.is_square(),
|
|
||||||
"Unable to compute the eigenvalue decomposition of a non-square matrix."
|
|
||||||
);
|
|
||||||
|
|
||||||
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 lda = n as i32;
|
|
||||||
|
|
||||||
// TODO: avoid the initialization?
|
|
||||||
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 info = 0;
|
|
||||||
let mut placeholder1 = [T::zero()];
|
|
||||||
let mut placeholder2 = [T::zero()];
|
|
||||||
|
|
||||||
let lwork = T::xgeev_work_size(
|
|
||||||
ljob,
|
ljob,
|
||||||
rjob,
|
rjob,
|
||||||
n as i32,
|
n as i32,
|
||||||
|
|
@ -274,142 +134,36 @@ where
|
||||||
lda,
|
lda,
|
||||||
wr.as_mut_slice(),
|
wr.as_mut_slice(),
|
||||||
wi.as_mut_slice(),
|
wi.as_mut_slice(),
|
||||||
&mut placeholder1,
|
vl_ref,
|
||||||
n as i32,
|
if left_eigenvectors { n as i32 } else { 1 },
|
||||||
&mut placeholder2,
|
vr_ref,
|
||||||
n as i32,
|
if eigenvectors { n as i32 } else { 1 },
|
||||||
|
&mut work,
|
||||||
|
lwork,
|
||||||
&mut info,
|
&mut info,
|
||||||
);
|
);
|
||||||
|
lapack_check!(info);
|
||||||
|
|
||||||
lapack_panic!(info);
|
Some(Self {
|
||||||
|
eigenvalues_re: wr,
|
||||||
|
eigenvalues_im: wi,
|
||||||
|
left_eigenvectors: vl,
|
||||||
|
eigenvectors: vr,
|
||||||
|
})
|
||||||
|
}
|
||||||
|
|
||||||
let mut work = vec![T::zero(); lwork as usize];
|
/// Returns `true` if all the eigenvalues are real.
|
||||||
|
pub fn eigenvalues_are_real(&self) -> bool {
|
||||||
match (left_eigenvectors, eigenvectors) {
|
self.eigenvalues_im.iter().all(|e| e.is_zero())
|
||||||
(true, true) => {
|
|
||||||
// TODO: avoid the initializations?
|
|
||||||
let mut vl = Matrix::zeros_generic(nrows, ncols);
|
|
||||||
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 vl.as_mut_slice(),
|
|
||||||
n 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());
|
|
||||||
}
|
|
||||||
(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());
|
|
||||||
}
|
|
||||||
(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());
|
|
||||||
}
|
|
||||||
(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());
|
|
||||||
}
|
|
||||||
(res, None, None)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/// The determinant of the decomposed matrix.
|
/// The determinant of the decomposed matrix.
|
||||||
#[inline]
|
#[inline]
|
||||||
#[must_use]
|
#[must_use]
|
||||||
pub fn determinant(&self) -> T {
|
pub fn determinant(&self) -> Complex<T> {
|
||||||
let mut det = T::one();
|
let mut det: Complex<T> = na::one();
|
||||||
for e in self.eigenvalues.iter() {
|
for (re, im) in self.eigenvalues_re.iter().zip(self.eigenvalues_im.iter()) {
|
||||||
det *= e.clone();
|
det *= Complex::new(re.clone(), im.clone());
|
||||||
}
|
}
|
||||||
|
|
||||||
det
|
det
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue