M-Labs/nalgebra
3
0
Fork 1
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Remove condition number, tests pass without. Add proper test generator for dynamic f64 type square matrices

Browse Source
This commit is contained in:
metric-space 2022-02-12 02:27:29 -05:00
parent 497a4d8bd9
commit fb0cb513e7
2 changed files with 58 additions and 104 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

98
nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs
View File

@ -1,74 +1,70 @@
use na::dimension::{Const, Dynamic};
use na::{DMatrix, EuclideanNorm, Norm, OMatrix};
use na::dimension::Const;
use na::{DMatrix, OMatrix};
use nl::GE;
use num_complex::Complex;
use simba::scalar::ComplexField;
use std::cmp;
use crate::proptest::*;
use proptest::{prop_assert, proptest};
use proptest::{prop_assert, prop_compose, proptest};
prop_compose! {
fn f64_squares() (n in PROPTEST_MATRIX_DIM) (a in matrix(PROPTEST_F64,n,n), b in matrix(PROPTEST_F64,n,n)) -> (DMatrix<f64>, DMatrix<f64>){
(a,b)
}}
proptest! {
#[test]
fn ge(n in PROPTEST_MATRIX_DIM) {
let n = cmp::max(1, cmp::min(n, 10));
let a = DMatrix::<f64>::new_random(n, n);
let b = DMatrix::<f64>::new_random(n, n);
let a_condition_no = a.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&a)));
let b_condition_no = b.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&b)));
fn ge((a,b) in f64_squares()){
if a_condition_no.unwrap_or(200000.0) < 5.0 && b_condition_no.unwrap_or(200000.0) < 5.0 {
let a_c = a.clone().map(|x| Complex::new(x, 0.0));
let b_c = b.clone().map(|x| Complex::new(x, 0.0));
let a_c = a.clone().map(|x| Complex::new(x, 0.0));
let b_c = b.clone().map(|x| Complex::new(x, 0.0));
let n = a.shape_generic().0;
let ge = GE::new(a.clone(), b.clone());
let (vsl,vsr) = ge.clone().eigenvectors();
let ge = GE::new(a.clone(), b.clone());
let (vsl,vsr) = ge.clone().eigenvectors();
for (i,(alpha,beta)) in ge.raw_eigenvalues().iter().enumerate() {
let l_a = a_c.clone() * Complex::new(*beta, 0.0);
let l_b = b_c.clone() * *alpha;
prop_assert!(
relative_eq!(
((&l_a - &l_b)*vsr.column(i)).map(|x| x.modulus()),
OMatrix::zeros_generic(Dynamic::new(n), Const::<1>),
epsilon = 1.0e-7));
for (i,(alpha,beta)) in ge.raw_eigenvalues().iter().enumerate() {
let l_a = a_c.clone() * Complex::new(*beta, 0.0);
let l_b = b_c.clone() * *alpha;
prop_assert!(
relative_eq!(
(vsl.column(i).adjoint()*(&l_a - &l_b)).map(|x| x.modulus()),
OMatrix::zeros_generic(Const::<1>, Dynamic::new(n)),
epsilon = 1.0e-7))
};
prop_assert!(
relative_eq!(
((&l_a - &l_b)*vsr.column(i)).map(|x| x.modulus()),
OMatrix::zeros_generic(n, Const::<1>),
epsilon = 1.0e-5));
prop_assert!(
relative_eq!(
(vsl.column(i).adjoint()*(&l_a - &l_b)).map(|x| x.modulus()),
OMatrix::zeros_generic(Const::<1>, n),
epsilon = 1.0e-5))
};
}
#[test]
fn ge_static(a in matrix4(), b in matrix4()) {
let a_condition_no = a.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&a)));
let b_condition_no = b.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&b)));
if a_condition_no.unwrap_or(200000.0) < 5.0 && b_condition_no.unwrap_or(200000.0) < 5.0 {
let ge = GE::new(a.clone(), b.clone());
let a_c =a.clone().map(|x| Complex::new(x, 0.0));
let b_c = b.clone().map(|x| Complex::new(x, 0.0));
let (vsl,vsr) = ge.eigenvectors();
let eigenvalues = ge.raw_eigenvalues();
let ge = GE::new(a.clone(), b.clone());
let a_c =a.clone().map(|x| Complex::new(x, 0.0));
let b_c = b.clone().map(|x| Complex::new(x, 0.0));
let (vsl,vsr) = ge.eigenvectors();
let eigenvalues = ge.raw_eigenvalues();
for (i,(alpha,beta)) in eigenvalues.iter().enumerate() {
let l_a = a_c.clone() * Complex::new(*beta, 0.0);
let l_b = b_c.clone() * *alpha;
for (i,(alpha,beta)) in eigenvalues.iter().enumerate() {
let l_a = a_c.clone() * Complex::new(*beta, 0.0);
let l_b = b_c.clone() * *alpha;
prop_assert!(
relative_eq!(
((&l_a - &l_b)*vsr.column(i)).map(|x| x.modulus()),
OMatrix::zeros_generic(Const::<4>, Const::<1>),
epsilon = 1.0e-7));
prop_assert!(
relative_eq!((vsl.column(i).adjoint()*(&l_a - &l_b)).map(|x| x.modulus()),
OMatrix::zeros_generic(Const::<1>, Const::<4>),
epsilon = 1.0e-7))
}
};
prop_assert!(
relative_eq!(
((&l_a - &l_b)*vsr.column(i)).map(|x| x.modulus()),
OMatrix::zeros_generic(Const::<4>, Const::<1>),
epsilon = 1.0e-5));
prop_assert!(
relative_eq!((vsl.column(i).adjoint()*(&l_a - &l_b)).map(|x| x.modulus()),
OMatrix::zeros_generic(Const::<1>, Const::<4>),
epsilon = 1.0e-5))
}
}
}

64
nalgebra-lapack/tests/linalg/qz.rs
View File

@ -1,74 +1,32 @@
use na::{DMatrix, EuclideanNorm, Norm};
use na::DMatrix;
use nl::QZ;
use num_complex::Complex;
use simba::scalar::ComplexField;
use std::cmp;
use crate::proptest::*;
use proptest::{prop_assert, proptest};
use proptest::{prop_assert, prop_compose, proptest};
prop_compose! {
fn f64_squares() (n in PROPTEST_MATRIX_DIM) (a in matrix(PROPTEST_F64,n,n), b in matrix(PROPTEST_F64,n,n)) -> (DMatrix<f64>, DMatrix<f64>){
(a,b)
}}
proptest! {
#[test]
fn qz(n in PROPTEST_MATRIX_DIM) {
let n = cmp::max(1, cmp::min(n, 10));
let a = DMatrix::<f64>::new_random(n, n);
let b = DMatrix::<f64>::new_random(n, n);
fn qz((a,b) in f64_squares()) {
let qz = QZ::new(a.clone(), b.clone());
let qz = QZ::new(a.clone(), b.clone());
let (vsl,s,t,vsr) = qz.clone().unpack();
let eigenvalues = qz.raw_eigenvalues();
prop_assert!(relative_eq!(&vsl * s * vsr.transpose(), a.clone(), epsilon = 1.0e-7));
prop_assert!(relative_eq!(vsl * t * vsr.transpose(), b.clone(), epsilon = 1.0e-7));
prop_assert!(relative_eq!(&vsl * s * vsr.transpose(), a, epsilon = 1.0e-7));
prop_assert!(relative_eq!(vsl * t * vsr.transpose(), b, epsilon = 1.0e-7));
let a_condition_no = a.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&a)));
let b_condition_no = b.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&b)));
if a_condition_no.unwrap_or(200000.0) < 5.0 && b_condition_no.unwrap_or(200000.0) < 5.0 {
let a_c = a.clone().map(|x| Complex::new(x, 0.0));
let b_c = b.clone().map(|x| Complex::new(x, 0.0));
for (alpha,beta) in eigenvalues.iter() {
let l_a = a_c.clone() * Complex::new(*beta, 0.0);
let l_b = b_c.clone() * *alpha;
prop_assert!(
relative_eq!(
(&l_a - &l_b).determinant().modulus(),
0.0,
epsilon = 1.0e-7));
};
};
}
#[test]
fn qz_static(a in matrix4(), b in matrix4()) {
let qz = QZ::new(a.clone(), b.clone());
let (vsl,s,t,vsr) = qz.unpack();
let eigenvalues = qz.raw_eigenvalues();
prop_assert!(relative_eq!(&vsl * s * vsr.transpose(), a, epsilon = 1.0e-7));
prop_assert!(relative_eq!(vsl * t * vsr.transpose(), b, epsilon = 1.0e-7));
let a_condition_no = a.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&a)));
let b_condition_no = b.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&b)));
if a_condition_no.unwrap_or(200000.0) < 5.0 && b_condition_no.unwrap_or(200000.0) < 5.0 {
let a_c =a.clone().map(|x| Complex::new(x, 0.0));
let b_c = b.clone().map(|x| Complex::new(x, 0.0));
for (alpha,beta) in eigenvalues.iter() {
let l_a = a_c.clone() * Complex::new(*beta, 0.0);
let l_b = b_c.clone() * *alpha;
prop_assert!(
relative_eq!(
(&l_a - &l_b).determinant().modulus(),
0.0,
epsilon = 1.0e-7));
}
};
}
}