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

cargo fmt + tests

Browse Source
This commit is contained in:
Sébastien Crozet 2022-10-30 17:22:08 +01:00
parent 9c8b5f0f38
commit e32f4ee16f
8 changed files with 184 additions and 112 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
nalgebra-lapack/examples/complex_eigen.rs
View File

@ -1,31 +0,0 @@
extern crate nalgebra as na;
extern crate nalgebra_lapack;
#[macro_use]
extern crate approx; // for assert_relative_eq
use na::Matrix3;
use nalgebra_lapack::Eigen;
use num_complex::Complex;
//Matrix taken from https://textbooks.math.gatech.edu/ila/1553/complex-eigenvalues.html
fn main() {
let m = Matrix3::<f64>::new(4.0/5.0, -3.0/5.0, 0.0, 3.0/5.0, 4.0/5.0, 0.0, 1.0, 2.0, 2.0);
let eigen = Eigen::new(m,true,true).expect("Eigen Creation Failed!");
let (some_eigenvalues, some_left_vec, some_right_vec) = eigen.get_complex_elements();
let eigenvalues = some_eigenvalues.expect("Eigenvalues Failed");
let _left_eigenvectors = some_left_vec.expect("Left Eigenvectors Failed");
let eigenvectors = some_right_vec.expect("Right Eigenvectors Failed");
assert_relative_eq!(eigenvalues[0].re, Complex::<f64>::new(4.0/5.0,3.0/5.0).re);
assert_relative_eq!(eigenvalues[0].im, Complex::<f64>::new(4.0/5.0,3.0/5.0).im);
assert_relative_eq!(eigenvalues[1].re, Complex::<f64>::new(4.0/5.0,-3.0/5.0).re);
assert_relative_eq!(eigenvalues[1].im, Complex::<f64>::new(4.0/5.0,-3.0/5.0).im);
assert_relative_eq!(eigenvectors[0][0].re, -12.0/32.7871926215100059134410999);
assert_relative_eq!(eigenvectors[0][0].im, -9.0/32.7871926215100059134410999);
assert_relative_eq!(eigenvectors[0][1].re, -9.0/32.7871926215100059134410999);
assert_relative_eq!(eigenvectors[0][1].im, 12.0/32.7871926215100059134410999);
assert_relative_eq!(eigenvectors[0][2].re, 25.0/32.7871926215100059134410999);
assert_relative_eq!(eigenvectors[0][2].im, 0.0);
}

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

@ -8,7 +8,7 @@ use simba::scalar::RealField;
use crate::ComplexHelper; use crate::ComplexHelper;
use na::dimension::{Const, Dim}; use na::dimension::{Const, Dim};
use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar, allocator::Allocator}; use na::{allocator::Allocator, DefaultAllocator, Matrix, OMatrix, OVector, Scalar};
use lapack; use lapack;
@ -147,7 +147,7 @@ where
eigenvalues_re: wr, eigenvalues_re: wr,
eigenvalues_im: wi, eigenvalues_im: wi,
left_eigenvectors: vl, left_eigenvectors: vl,
eigenvectors: vr eigenvectors: vr,
}) })
} }
@ -169,17 +169,30 @@ where
} }
/// Returns a tuple of vectors. The elements of the tuple are the real parts of the eigenvalues, left eigenvectors and right eigenvectors respectively. /// Returns a tuple of vectors. The elements of the tuple are the real parts of the eigenvalues, left eigenvectors and right eigenvectors respectively.
pub fn get_real_elements(&self) -> (Vec<T>, Option<Vec<OVector<T, D>>>, Option<Vec<OVector<T, D>>>) where DefaultAllocator: Allocator<T, D> { pub fn get_real_elements(
&self,
) -> (
Vec<T>,
Option<Vec<OVector<T, D>>>,
Option<Vec<OVector<T, D>>>,
)
where
DefaultAllocator: Allocator<T, D>,
{
let (number_of_elements, _) = self.eigenvalues_re.shape_generic(); let (number_of_elements, _) = self.eigenvalues_re.shape_generic();
let number_of_elements_value = number_of_elements.value(); let number_of_elements_value = number_of_elements.value();
let mut eigenvalues = Vec::<T>::with_capacity(number_of_elements_value); let mut eigenvalues = Vec::<T>::with_capacity(number_of_elements_value);
let mut eigenvectors = match self.eigenvectors.is_some() { let mut eigenvectors = match self.eigenvectors.is_some() {
true => Some(Vec::<OVector<T, D>>::with_capacity(number_of_elements_value)), true => Some(Vec::<OVector<T, D>>::with_capacity(
false => None number_of_elements_value,
)),
false => None,
}; };
let mut left_eigenvectors = match self.left_eigenvectors.is_some() { let mut left_eigenvectors = match self.left_eigenvectors.is_some() {
true => Some(Vec::<OVector<T, D>>::with_capacity(number_of_elements_value)), true => Some(Vec::<OVector<T, D>>::with_capacity(
false => None number_of_elements_value,
)),
false => None,
}; };
let mut c = 0; let mut c = 0;
@ -187,11 +200,21 @@ where
eigenvalues.push(self.eigenvalues_re[c].clone()); eigenvalues.push(self.eigenvalues_re[c].clone());
if eigenvectors.is_some() { if eigenvectors.is_some() {
eigenvectors.as_mut().unwrap().push((&self.eigenvectors.as_ref()).unwrap().column(c).into_owned()); eigenvectors.as_mut().unwrap().push(
(&self.eigenvectors.as_ref())
.unwrap()
.column(c)
.into_owned(),
);
} }
if left_eigenvectors.is_some() { if left_eigenvectors.is_some() {
left_eigenvectors.as_mut().unwrap().push((&self.left_eigenvectors.as_ref()).unwrap().column(c).into_owned()); left_eigenvectors.as_mut().unwrap().push(
(&self.left_eigenvectors.as_ref())
.unwrap()
.column(c)
.into_owned(),
);
} }
if self.eigenvalues_im[c] != T::zero() { if self.eigenvalues_im[c] != T::zero() {
//skip next entry //skip next entry
@ -200,43 +223,75 @@ where
c += 1; c += 1;
} }
(eigenvalues, left_eigenvectors, eigenvectors) (eigenvalues, left_eigenvectors, eigenvectors)
} }
/// Returns a tuple of vectors. The elements of the tuple are the complex eigenvalues, complex left eigenvectors and complex right eigenvectors respectively. /// Returns a tuple of vectors. The elements of the tuple are the complex eigenvalues, complex left eigenvectors and complex right eigenvectors respectively.
/// The elements appear as conjugate pairs within each vector, with the positive of the pair always being first. /// The elements appear as conjugate pairs within each vector, with the positive of the pair always being first.
pub fn get_complex_elements(&self) -> (Option<Vec<Complex<T>>>, Option<Vec<OVector<Complex<T>, D>>>, Option<Vec<OVector<Complex<T>, D>>>) where DefaultAllocator: Allocator<Complex<T>, D> { pub fn get_complex_elements(
&self,
) -> (
Option<Vec<Complex<T>>>,
Option<Vec<OVector<Complex<T>, D>>>,
Option<Vec<OVector<Complex<T>, D>>>,
)
where
DefaultAllocator: Allocator<Complex<T>, D>,
{
match self.eigenvalues_are_real() { match self.eigenvalues_are_real() {
true => (None, None, None), true => (None, None, None),
false => { false => {
let (number_of_elements, _) = self.eigenvalues_re.shape_generic(); let (number_of_elements, _) = self.eigenvalues_re.shape_generic();
let number_of_elements_value = number_of_elements.value(); let number_of_elements_value = number_of_elements.value();
let number_of_complex_entries = self.eigenvalues_im.iter().fold(0, |acc, e| if !e.is_zero() {acc + 1} else {acc}); let number_of_complex_entries =
self.eigenvalues_im
.iter()
.fold(0, |acc, e| if !e.is_zero() { acc + 1 } else { acc });
let mut eigenvalues = Vec::<Complex<T>>::with_capacity(number_of_complex_entries); let mut eigenvalues = Vec::<Complex<T>>::with_capacity(number_of_complex_entries);
let mut eigenvectors = match self.eigenvectors.is_some() { let mut eigenvectors = match self.eigenvectors.is_some() {
true => Some(Vec::<OVector<Complex<T>, D>>::with_capacity(number_of_complex_entries)), true => Some(Vec::<OVector<Complex<T>, D>>::with_capacity(
false => None number_of_complex_entries,
)),
false => None,
}; };
let mut left_eigenvectors = match self.left_eigenvectors.is_some() { let mut left_eigenvectors = match self.left_eigenvectors.is_some() {
true => Some(Vec::<OVector<Complex<T>, D>>::with_capacity(number_of_complex_entries)), true => Some(Vec::<OVector<Complex<T>, D>>::with_capacity(
false => None number_of_complex_entries,
)),
false => None,
}; };
let mut c = 0; let mut c = 0;
while c < number_of_elements_value { while c < number_of_elements_value {
if self.eigenvalues_im[c] != T::zero() { if self.eigenvalues_im[c] != T::zero() {
//Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first. //Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first.
eigenvalues.push(Complex::<T>::new(self.eigenvalues_re[c].clone(), self.eigenvalues_im[c].clone())); eigenvalues.push(Complex::<T>::new(
eigenvalues.push(Complex::<T>::new(self.eigenvalues_re[c+1].clone(), self.eigenvalues_im[c+1].clone())); self.eigenvalues_re[c].clone(),
self.eigenvalues_im[c].clone(),
));
eigenvalues.push(Complex::<T>::new(
self.eigenvalues_re[c + 1].clone(),
self.eigenvalues_im[c + 1].clone(),
));
if eigenvectors.is_some() { if eigenvectors.is_some() {
let mut vec = OVector::<Complex<T>, D>::zeros_generic(number_of_elements, Const::<1>); let mut vec = OVector::<Complex<T>, D>::zeros_generic(
let mut vec_conj = OVector::<Complex<T>, D>::zeros_generic(number_of_elements, Const::<1>); number_of_elements,
Const::<1>,
);
let mut vec_conj = OVector::<Complex<T>, D>::zeros_generic(
number_of_elements,
Const::<1>,
);
for r in 0..number_of_elements_value { for r in 0..number_of_elements_value {
vec[r] = Complex::<T>::new((&self.eigenvectors.as_ref()).unwrap()[(r,c)].clone(),(&self.eigenvectors.as_ref()).unwrap()[(r,c+1)].clone()); vec[r] = Complex::<T>::new(
vec_conj[r] = Complex::<T>::new((&self.eigenvectors.as_ref()).unwrap()[(r,c)].clone(),(&self.eigenvectors.as_ref()).unwrap()[(r,c+1)].clone()); (&self.eigenvectors.as_ref()).unwrap()[(r, c)].clone(),
(&self.eigenvectors.as_ref()).unwrap()[(r, c + 1)].clone(),
);
vec_conj[r] = Complex::<T>::new(
(&self.eigenvectors.as_ref()).unwrap()[(r, c)].clone(),
(&self.eigenvectors.as_ref()).unwrap()[(r, c + 1)].clone(),
);
} }
eigenvectors.as_mut().unwrap().push(vec); eigenvectors.as_mut().unwrap().push(vec);
@ -244,12 +299,24 @@ where
} }
if left_eigenvectors.is_some() { if left_eigenvectors.is_some() {
let mut vec = OVector::<Complex<T>, D>::zeros_generic(number_of_elements, Const::<1>); let mut vec = OVector::<Complex<T>, D>::zeros_generic(
let mut vec_conj = OVector::<Complex<T>, D>::zeros_generic(number_of_elements, Const::<1>); number_of_elements,
Const::<1>,
);
let mut vec_conj = OVector::<Complex<T>, D>::zeros_generic(
number_of_elements,
Const::<1>,
);
for r in 0..number_of_elements_value { for r in 0..number_of_elements_value {
vec[r] = Complex::<T>::new((&self.left_eigenvectors.as_ref()).unwrap()[(r,c)].clone(),(&self.left_eigenvectors.as_ref()).unwrap()[(r,c+1)].clone()); vec[r] = Complex::<T>::new(
vec_conj[r] = Complex::<T>::new((&self.left_eigenvectors.as_ref()).unwrap()[(r,c)].clone(),(&self.left_eigenvectors.as_ref()).unwrap()[(r,c+1)].clone()); (&self.left_eigenvectors.as_ref()).unwrap()[(r, c)].clone(),
(&self.left_eigenvectors.as_ref()).unwrap()[(r, c + 1)].clone(),
);
vec_conj[r] = Complex::<T>::new(
(&self.left_eigenvectors.as_ref()).unwrap()[(r, c)].clone(),
(&self.left_eigenvectors.as_ref()).unwrap()[(r, c + 1)].clone(),
);
} }
left_eigenvectors.as_mut().unwrap().push(vec); left_eigenvectors.as_mut().unwrap().push(vec);
@ -264,7 +331,6 @@ where
} }
} }
} }
} }
/* /*

6
nalgebra-lapack/tests/linalg/cholesky.rs
View File

@ -58,8 +58,8 @@ proptest! {
let sol1 = chol.solve(&b1).unwrap(); let sol1 = chol.solve(&b1).unwrap();
let sol2 = chol.solve(&b2).unwrap(); let sol2 = chol.solve(&b2).unwrap();
prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-7)); prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-4));
prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-7)); prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-4));
} }
} }
@ -84,7 +84,7 @@ proptest! {
let id1 = &m * &m1; let id1 = &m * &m1;
let id2 = &m1 * &m; let id2 = &m1 * &m;
prop_assert!(id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)) prop_assert!(id1.is_identity(1.0e-4) && id2.is_identity(1.0e-4))
} }
} }
} }

56
nalgebra-lapack/tests/linalg/complex_eigen.rs
View File

@ -1,19 +1,47 @@
use std::cmp; use na::Matrix3;
use na::{Matrix3};
use nalgebra_lapack::Eigen; use nalgebra_lapack::Eigen;
use num_complex::Complex;
use crate::proptest::*; #[test]
use proptest::{prop_assert, proptest}; fn complex_eigen() {
let m = Matrix3::<f64>::new(
4.0 / 5.0,
-3.0 / 5.0,
0.0,
3.0 / 5.0,
4.0 / 5.0,
0.0,
1.0,
2.0,
2.0,
);
let eigen = Eigen::new(m, true, true).expect("Eigen Creation Failed!");
let (some_eigenvalues, some_left_vec, some_right_vec) = eigen.get_complex_elements();
let eigenvalues = some_eigenvalues.expect("Eigenvalues Failed");
let _left_eigenvectors = some_left_vec.expect("Left Eigenvectors Failed");
let eigenvectors = some_right_vec.expect("Right Eigenvectors Failed");
proptest! { assert_relative_eq!(
//#[test] eigenvalues[0].re,
// fn complex_eigen() { Complex::<f64>::new(4.0 / 5.0, 3.0 / 5.0).re
// let n = cmp::max(1, cmp::min(n, 10)); );
// let m = DMatrix::<f64>::new_random(n, n); assert_relative_eq!(
// let eig = SymmetricEigen::new(m.clone()); eigenvalues[0].im,
// let recomp = eig.recompose(); Complex::<f64>::new(4.0 / 5.0, 3.0 / 5.0).im
// prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)) );
// } assert_relative_eq!(
eigenvalues[1].re,
Complex::<f64>::new(4.0 / 5.0, -3.0 / 5.0).re
);
assert_relative_eq!(
eigenvalues[1].im,
Complex::<f64>::new(4.0 / 5.0, -3.0 / 5.0).im
);
assert_relative_eq!(eigenvectors[0][0].re, -12.0 / 32.7871926215100059134410999);
assert_relative_eq!(eigenvectors[0][0].im, -9.0 / 32.7871926215100059134410999);
assert_relative_eq!(eigenvectors[0][1].re, -9.0 / 32.7871926215100059134410999);
assert_relative_eq!(eigenvectors[0][1].im, 12.0 / 32.7871926215100059134410999);
assert_relative_eq!(eigenvectors[0][2].re, 25.0 / 32.7871926215100059134410999);
assert_relative_eq!(eigenvectors[0][2].im, 0.0);
} }

16
nalgebra-lapack/tests/linalg/lu.rs
View File

@ -51,10 +51,10 @@ proptest! {
let tr_sol1 = lup.solve_transpose(&b1).unwrap(); let tr_sol1 = lup.solve_transpose(&b1).unwrap();
let tr_sol2 = lup.solve_transpose(&b2).unwrap(); let tr_sol2 = lup.solve_transpose(&b2).unwrap();
prop_assert!(relative_eq!(&m * sol1, b1, epsilon = 1.0e-7)); prop_assert!(relative_eq!(&m * sol1, b1, epsilon = 1.0e-5));
prop_assert!(relative_eq!(&m * sol2, b2, epsilon = 1.0e-7)); prop_assert!(relative_eq!(&m * sol2, b2, epsilon = 1.0e-5));
prop_assert!(relative_eq!(m.transpose() * tr_sol1, b1, epsilon = 1.0e-7)); prop_assert!(relative_eq!(m.transpose() * tr_sol1, b1, epsilon = 1.0e-5));
prop_assert!(relative_eq!(m.transpose() * tr_sol2, b2, epsilon = 1.0e-7)); prop_assert!(relative_eq!(m.transpose() * tr_sol2, b2, epsilon = 1.0e-5));
} }
#[test] #[test]
@ -68,10 +68,10 @@ proptest! {
let tr_sol1 = lup.solve_transpose(&b1).unwrap(); let tr_sol1 = lup.solve_transpose(&b1).unwrap();
let tr_sol2 = lup.solve_transpose(&b2).unwrap(); let tr_sol2 = lup.solve_transpose(&b2).unwrap();
prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-7)); prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-5));
prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-7)); prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-5));
prop_assert!(relative_eq!(m.transpose() * tr_sol1, b1, epsilon = 1.0e-7)); prop_assert!(relative_eq!(m.transpose() * tr_sol1, b1, epsilon = 1.0e-5));
prop_assert!(relative_eq!(m.transpose() * tr_sol2, b2, epsilon = 1.0e-7)); prop_assert!(relative_eq!(m.transpose() * tr_sol2, b2, epsilon = 1.0e-5));
} }
#[test] #[test]

2
nalgebra-lapack/tests/linalg/mod.rs
View File

@ -1,4 +1,5 @@
mod cholesky; mod cholesky;
mod complex_eigen;
mod generalized_eigenvalues; mod generalized_eigenvalues;
mod lu; mod lu;
mod qr; mod qr;
@ -7,4 +8,3 @@ mod real_eigensystem;
mod schur; mod schur;
mod svd; mod svd;
mod symmetric_eigen; mod symmetric_eigen;
mod complex_eigen;

18
nalgebra-lapack/tests/linalg/real_eigensystem.rs
View File

@ -13,30 +13,36 @@ proptest! {
let m = DMatrix::<f64>::new_random(n, n); let m = DMatrix::<f64>::new_random(n, n);
if let Some(eig) = Eigen::new(m.clone(), true, true) { if let Some(eig) = Eigen::new(m.clone(), true, true) {
let eigvals = DMatrix::from_diagonal(&eig.eigenvalues); // TODO: test the complex case too.
if eig.eigenvalues_are_real() {
let eigvals = DMatrix::from_diagonal(&eig.eigenvalues_re);
let transformed_eigvectors = &m * eig.eigenvectors.as_ref().unwrap(); let transformed_eigvectors = &m * eig.eigenvectors.as_ref().unwrap();
let scaled_eigvectors = eig.eigenvectors.as_ref().unwrap() * &eigvals; let scaled_eigvectors = eig.eigenvectors.as_ref().unwrap() * &eigvals;
let transformed_left_eigvectors = m.transpose() * eig.left_eigenvectors.as_ref().unwrap(); let transformed_left_eigvectors = m.transpose() * eig.left_eigenvectors.as_ref().unwrap();
let scaled_left_eigvectors = eig.left_eigenvectors.as_ref().unwrap() * &eigvals; let scaled_left_eigvectors = eig.left_eigenvectors.as_ref().unwrap() * &eigvals;
prop_assert!(relative_eq!(transformed_eigvectors, scaled_eigvectors, epsilon = 1.0e-7)); prop_assert!(relative_eq!(transformed_eigvectors, scaled_eigvectors, epsilon = 1.0e-5));
prop_assert!(relative_eq!(transformed_left_eigvectors, scaled_left_eigvectors, epsilon = 1.0e-7)); prop_assert!(relative_eq!(transformed_left_eigvectors, scaled_left_eigvectors, epsilon = 1.0e-5));
}
} }
} }
#[test] #[test]
fn eigensystem_static(m in matrix4()) { fn eigensystem_static(m in matrix4()) {
if let Some(eig) = Eigen::new(m, true, true) { if let Some(eig) = Eigen::new(m, true, true) {
let eigvals = Matrix4::from_diagonal(&eig.eigenvalues); // TODO: test the complex case too.
if eig.eigenvalues_are_real() {
let eigvals = Matrix4::from_diagonal(&eig.eigenvalues_re);
let transformed_eigvectors = m * eig.eigenvectors.unwrap(); let transformed_eigvectors = m * eig.eigenvectors.unwrap();
let scaled_eigvectors = eig.eigenvectors.unwrap() * eigvals; let scaled_eigvectors = eig.eigenvectors.unwrap() * eigvals;
let transformed_left_eigvectors = m.transpose() * eig.left_eigenvectors.unwrap(); let transformed_left_eigvectors = m.transpose() * eig.left_eigenvectors.unwrap();
let scaled_left_eigvectors = eig.left_eigenvectors.unwrap() * eigvals; let scaled_left_eigvectors = eig.left_eigenvectors.unwrap() * eigvals;
prop_assert!(relative_eq!(transformed_eigvectors, scaled_eigvectors, epsilon = 1.0e-7)); prop_assert!(relative_eq!(transformed_eigvectors, scaled_eigvectors, epsilon = 1.0e-5));
prop_assert!(relative_eq!(transformed_left_eigvectors, scaled_left_eigvectors, epsilon = 1.0e-7)); prop_assert!(relative_eq!(transformed_left_eigvectors, scaled_left_eigvectors, epsilon = 1.0e-5));
}
} }
} }
} }

13
nalgebra-lapack/tests/linalg/schur.rs
View File

@ -11,14 +11,17 @@ proptest! {
let n = cmp::max(1, cmp::min(n, 10)); let n = cmp::max(1, cmp::min(n, 10));
let m = DMatrix::<f64>::new_random(n, n); let m = DMatrix::<f64>::new_random(n, n);
let (vecs, vals) = Schur::new(m.clone()).unpack(); if let Some(schur) = Schur::try_new(m.clone()) {
let (vecs, vals) = schur.unpack();
prop_assert!(relative_eq!(&vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7)) prop_assert!(relative_eq!(&vecs * vals * vecs.transpose(), m, epsilon = 1.0e-5))
}
} }
#[test] #[test]
fn schur_static(m in matrix4()) { fn schur_static(m in matrix4()) {
let (vecs, vals) = Schur::new(m.clone()).unpack(); if let Some(schur) = Schur::try_new(m.clone()) {
prop_assert!(relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7)) let (vecs, vals) = schur.unpack();
prop_assert!(relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-5))
}
} }
} }