forked from M-Labs/nalgebra
180 lines
8.7 KiB
Rust
180 lines
8.7 KiB
Rust
use std::cmp;
|
|
|
|
use na::{DMatrix, Matrix2, Matrix3, Matrix4, SymmetricEigen};
|
|
|
|
#[cfg(feature = "arbitrary")]
|
|
quickcheck! {
|
|
fn symmetric_eigen(n: usize) -> bool {
|
|
let n = cmp::max(1, cmp::min(n, 10));
|
|
let m = DMatrix::<f64>::new_random(n, n);
|
|
let eig = SymmetricEigen::new(m.clone());
|
|
let recomp = eig.recompose();
|
|
|
|
println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
|
|
|
|
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
|
|
}
|
|
|
|
fn symmetric_eigen_singular(n: usize) -> bool {
|
|
let n = cmp::max(1, cmp::min(n, 10));
|
|
let mut m = DMatrix::<f64>::new_random(n, n);
|
|
m.row_mut(n / 2).fill(0.0);
|
|
m.column_mut(n / 2).fill(0.0);
|
|
let eig = SymmetricEigen::new(m.clone());
|
|
let recomp = eig.recompose();
|
|
|
|
println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
|
|
|
|
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
|
|
}
|
|
|
|
fn symmetric_eigen_static_square_4x4(m: Matrix4<f64>) -> bool {
|
|
let eig = SymmetricEigen::new(m);
|
|
let recomp = eig.recompose();
|
|
|
|
println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
|
|
|
|
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
|
|
}
|
|
|
|
fn symmetric_eigen_static_square_3x3(m: Matrix3<f64>) -> bool {
|
|
let eig = SymmetricEigen::new(m);
|
|
let recomp = eig.recompose();
|
|
|
|
println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
|
|
|
|
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
|
|
}
|
|
|
|
fn symmetric_eigen_static_square_2x2(m: Matrix2<f64>) -> bool {
|
|
let eig = SymmetricEigen::new(m);
|
|
let recomp = eig.recompose();
|
|
|
|
println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
|
|
|
|
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
|
|
}
|
|
}
|
|
|
|
// Test proposed on the issue #176 of rulinalg.
|
|
#[test]
|
|
fn symmetric_eigen_singular_24x24() {
|
|
let m = DMatrix::from_row_slice(24, 24, &[
|
|
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
-1.0, -1.0, -1.0, -1.0, -1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 0.0, 1.0, 1.0, 1.0,
|
|
0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0]);
|
|
|
|
let eig = SymmetricEigen::new(m.clone());
|
|
let recomp = eig.recompose();
|
|
|
|
assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5));
|
|
}
|
|
|
|
|
|
// #[cfg(feature = "arbitrary")]
|
|
// quickcheck! {
|
|
// FIXME: full eigendecomposition is not implemented yet because of its complexity when some
|
|
// eigenvalues have multiplicity > 1.
|
|
//
|
|
// /*
|
|
// * NOTE: for the following tests, we use only upper-triangular matrices.
|
|
// * Thes ensures the schur decomposition will work, and allows use to test the eigenvector
|
|
// * computation.
|
|
// */
|
|
// fn eigen(n: usize) -> bool {
|
|
// let n = cmp::max(1, cmp::min(n, 10));
|
|
// let m = DMatrix::<f64>::new_random(n, n).upper_triangle();
|
|
//
|
|
// let eig = RealEigen::new(m.clone()).unwrap();
|
|
// verify_eigenvectors(m, eig)
|
|
// }
|
|
//
|
|
// fn eigen_with_adjascent_duplicate_diagonals(n: usize) -> bool {
|
|
// let n = cmp::max(1, cmp::min(n, 10));
|
|
// let mut m = DMatrix::<f64>::new_random(n, n).upper_triangle();
|
|
//
|
|
// // Suplicate some adjascent diagonal elements.
|
|
// for i in 0 .. n / 2 {
|
|
// m[(i * 2 + 1, i * 2 + 1)] = m[(i * 2, i * 2)];
|
|
// }
|
|
//
|
|
// let eig = RealEigen::new(m.clone()).unwrap();
|
|
// verify_eigenvectors(m, eig)
|
|
// }
|
|
//
|
|
// fn eigen_with_nonadjascent_duplicate_diagonals(n: usize) -> bool {
|
|
// let n = cmp::max(3, cmp::min(n, 10));
|
|
// let mut m = DMatrix::<f64>::new_random(n, n).upper_triangle();
|
|
//
|
|
// // Suplicate some diagonal elements.
|
|
// for i in n / 2 .. n {
|
|
// m[(i, i)] = m[(i - n / 2, i - n / 2)];
|
|
// }
|
|
//
|
|
// let eig = RealEigen::new(m.clone()).unwrap();
|
|
// verify_eigenvectors(m, eig)
|
|
// }
|
|
//
|
|
// fn eigen_static_square_4x4(m: Matrix4<f64>) -> bool {
|
|
// let m = m.upper_triangle();
|
|
// let eig = RealEigen::new(m.clone()).unwrap();
|
|
// verify_eigenvectors(m, eig)
|
|
// }
|
|
//
|
|
// fn eigen_static_square_3x3(m: Matrix3<f64>) -> bool {
|
|
// let m = m.upper_triangle();
|
|
// let eig = RealEigen::new(m.clone()).unwrap();
|
|
// verify_eigenvectors(m, eig)
|
|
// }
|
|
//
|
|
// fn eigen_static_square_2x2(m: Matrix2<f64>) -> bool {
|
|
// let m = m.upper_triangle();
|
|
// println!("{}", m);
|
|
// let eig = RealEigen::new(m.clone()).unwrap();
|
|
// verify_eigenvectors(m, eig)
|
|
// }
|
|
// }
|
|
//
|
|
// fn verify_eigenvectors<D: Dim>(m: MatrixN<f64, D>, mut eig: RealEigen<f64, D>) -> bool
|
|
// where DefaultAllocator: Allocator<f64, D, D> +
|
|
// Allocator<f64, D> +
|
|
// Allocator<usize, D, D> +
|
|
// Allocator<usize, D>,
|
|
// MatrixN<f64, D>: Display,
|
|
// VectorN<f64, D>: Display {
|
|
// let mv = &m * &eig.eigenvectors;
|
|
//
|
|
// println!("eigenvalues: {}eigenvectors: {}", eig.eigenvalues, eig.eigenvectors);
|
|
//
|
|
// let dim = m.nrows();
|
|
// for i in 0 .. dim {
|
|
// let mut col = eig.eigenvectors.column_mut(i);
|
|
// col *= eig.eigenvalues[i];
|
|
// }
|
|
//
|
|
// println!("{}{:.5}{:.5}", m, mv, eig.eigenvectors);
|
|
//
|
|
// relative_eq!(eig.eigenvectors, mv, epsilon = 1.0e-5)
|
|
// }
|