95b021dbfb
Don't enable `serde/std` by default. Add a `serde-serialize-no-std` feature to enable serde without its default features.
115 lines
3.4 KiB
Rust
115 lines
3.4 KiB
Rust
#[cfg(feature = "serde-serialize-no-std")]
|
|
use serde::{Deserialize, Serialize};
|
|
|
|
use num_complex::Complex;
|
|
use simba::scalar::ComplexField;
|
|
use std::cmp;
|
|
use std::fmt::Display;
|
|
use std::ops::MulAssign;
|
|
|
|
use crate::allocator::Allocator;
|
|
use crate::base::dimension::{Dim, DimDiff, DimSub, Dynamic, U1, U2, U3};
|
|
use crate::base::storage::Storage;
|
|
use crate::base::{
|
|
DefaultAllocator, Hessenberg, OMatrix, OVector, SquareMatrix, Unit, Vector2, Vector3,
|
|
};
|
|
use crate::constraint::{DimEq, ShapeConstraint};
|
|
|
|
use crate::geometry::{Reflection, UnitComplex};
|
|
use crate::linalg::householder;
|
|
use crate::linalg::Schur;
|
|
|
|
/// Eigendecomposition of a real matrix with real eigenvalues (or complex eigen values for complex matrices).
|
|
#[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))]
|
|
#[cfg_attr(
|
|
feature = "serde-serialize-no-std",
|
|
serde(bound(serialize = "DefaultAllocator: Allocator<T, D>,
|
|
OVector<T, D>: Serialize,
|
|
OMatrix<T, D, D>: Serialize"))
|
|
)]
|
|
#[cfg_attr(
|
|
feature = "serde-serialize-no-std",
|
|
serde(bound(deserialize = "DefaultAllocator: Allocator<T, D>,
|
|
OVector<T, D>: Serialize,
|
|
OMatrix<T, D, D>: Deserialize<'de>"))
|
|
)]
|
|
#[derive(Clone, Debug)]
|
|
pub struct Eigen<T: ComplexField, D: Dim>
|
|
where
|
|
DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
|
|
{
|
|
pub eigenvectors: OMatrix<T, D, D>,
|
|
pub eigenvalues: OVector<T, D>,
|
|
}
|
|
|
|
impl<T: ComplexField, D: Dim> Copy for Eigen<T, D>
|
|
where
|
|
DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
|
|
OMatrix<T, D, D>: Copy,
|
|
OVector<T, D>: Copy,
|
|
{
|
|
}
|
|
|
|
impl<T: ComplexField, D: Dim> Eigen<T, D>
|
|
where
|
|
D: DimSub<U1>, // For Hessenberg.
|
|
ShapeConstraint: DimEq<Dynamic, DimDiff<D, U1>>, // For Hessenberg.
|
|
DefaultAllocator: Allocator<T, D, DimDiff<D, U1>>
|
|
+ Allocator<T, DimDiff<D, U1>>
|
|
+ Allocator<T, D, D>
|
|
+ Allocator<T, D>,
|
|
// XXX: for debug
|
|
DefaultAllocator: Allocator<usize, D, D>,
|
|
OMatrix<T, D, D>: Display,
|
|
{
|
|
/// Computes the eigendecomposition of a diagonalizable matrix with Complex eigenvalues.
|
|
pub fn new(m: OMatrix<T, D, D>) -> Option<Eigen<T, D>> {
|
|
assert!(
|
|
m.is_square(),
|
|
"Unable to compute the eigendecomposition of a non-square matrix."
|
|
);
|
|
|
|
let dim = m.nrows();
|
|
let (mut eigenvectors, mut eigenvalues) = Schur::new(m, 0).unwrap().unpack();
|
|
|
|
println!("Schur eigenvalues: {}", eigenvalues);
|
|
|
|
// Check that the eigenvalues are all Complex.
|
|
for i in 0..dim - 1 {
|
|
if !eigenvalues[(i + 1, i)].is_zero() {
|
|
return None;
|
|
}
|
|
}
|
|
|
|
for j in 1..dim {
|
|
for i in 0..j {
|
|
let diff = eigenvalues[(i, i)] - eigenvalues[(j, j)];
|
|
|
|
if diff.is_zero() && !eigenvalues[(i, j)].is_zero() {
|
|
return None;
|
|
}
|
|
|
|
let z = -eigenvalues[(i, j)] / diff;
|
|
|
|
for k in j + 1..dim {
|
|
eigenvalues[(i, k)] -= z * eigenvalues[(j, k)];
|
|
}
|
|
|
|
for k in 0..dim {
|
|
eigenvectors[(k, j)] += z * eigenvectors[(k, i)];
|
|
}
|
|
}
|
|
}
|
|
|
|
// Normalize the eigenvector basis.
|
|
for i in 0..dim {
|
|
let _ = eigenvectors.column_mut(i).normalize_mut();
|
|
}
|
|
|
|
Some(Eigen {
|
|
eigenvectors,
|
|
eigenvalues: eigenvalues.diagonal(),
|
|
})
|
|
}
|
|
}
|