Restructure test modules to avoid warnings

These warnings occurred only when running the test suite with no
features. Lots of uses had to be rescoped into newly created modules to
make it easier to separate these issues.
This commit is contained in:
Eduard Bopp 2018-01-17 16:48:47 +01:00
parent 4a926736fe
commit bba1993e58
17 changed files with 673 additions and 649 deletions

View File

@ -1,9 +1,8 @@
use std::cmp; #![cfg(feature = "arbitrary")]
use std::cmp;
use na::{DVector, DMatrix}; use na::{DVector, DMatrix};
#[cfg(feature = "arbitrary")]
quickcheck! { quickcheck! {
/* /*
* *

View File

@ -1,8 +1,5 @@
use num::{Zero, One}; use num::{Zero, One};
use num::Float; use num::Float;
use std::fmt::Display;
use alga::linear::FiniteDimInnerSpace;
use na::{self, use na::{self,
DVector, DMatrix, DVector, DMatrix,
@ -806,93 +803,100 @@ mod normalization_tests {
} }
} }
#[cfg(feature = "arbitrary")]
// FIXME: move this to alga ? // FIXME: move this to alga ?
macro_rules! finite_dim_inner_space_test( mod finite_dim_inner_space_tests {
($($Vector: ident, $orthonormal_subspace: ident, $orthonormalization: ident);* $(;)*) => {$( use super::*;
#[cfg(feature = "arbitrary")] use std::fmt::Display;
quickcheck!{ use alga::linear::FiniteDimInnerSpace;
fn $orthonormal_subspace(vs: Vec<$Vector<f64>>) -> bool {
let mut given_basis = vs.clone();
let given_basis_dim = $Vector::orthonormalize(&mut given_basis[..]);
let mut ortho_basis = Vec::new();
$Vector::orthonormal_subspace_basis(
&given_basis[.. given_basis_dim],
|e| { ortho_basis.push(*e); true }
);
if !is_subspace_basis(&ortho_basis[..]) { macro_rules! finite_dim_inner_space_test(
return false; ($($Vector: ident, $orthonormal_subspace: ident, $orthonormalization: ident);* $(;)*) => {$(
quickcheck!{
fn $orthonormal_subspace(vs: Vec<$Vector<f64>>) -> bool {
let mut given_basis = vs.clone();
let given_basis_dim = $Vector::orthonormalize(&mut given_basis[..]);
let mut ortho_basis = Vec::new();
$Vector::orthonormal_subspace_basis(
&given_basis[.. given_basis_dim],
|e| { ortho_basis.push(*e); true }
);
if !is_subspace_basis(&ortho_basis[..]) {
return false;
}
for v in vs {
for b in &ortho_basis {
if !relative_eq!(v.dot(b), 0.0, epsilon = 1.0e-7) {
println!("Found dot product: {} · {} = {}", v, b, v.dot(b));
return false;
}
}
}
true
} }
for v in vs { fn $orthonormalization(vs: Vec<$Vector<f64>>) -> bool {
for b in &ortho_basis { let mut basis = vs.clone();
if !relative_eq!(v.dot(b), 0.0, epsilon = 1.0e-7) { let subdim = $Vector::orthonormalize(&mut basis[..]);
println!("Found dot product: {} · {} = {}", v, b, v.dot(b));
if !is_subspace_basis(&basis[.. subdim]) {
return false;
}
for mut e in vs {
for b in &basis[.. subdim] {
e -= e.dot(b) * b
}
// Any element of `e` must be a linear combination of the basis elements.
if !relative_eq!(e.norm(), 0.0, epsilon = 1.0e-7) {
println!("Orthonormalization; element decomposition failure: {}", e);
println!("... the non-zero norm is: {}", e.norm());
return false; return false;
} }
} }
}
true true
}
}
)*}
);
finite_dim_inner_space_test!(
Vector1, orthonormal_subspace_basis1, orthonormalize1;
Vector2, orthonormal_subspace_basis2, orthonormalize2;
Vector3, orthonormal_subspace_basis3, orthonormalize3;
Vector4, orthonormal_subspace_basis4, orthonormalize4;
Vector5, orthonormal_subspace_basis5, orthonormalize5;
Vector6, orthonormal_subspace_basis6, orthonormalize6;
);
/*
*
* Helper functions.
*
*/
#[cfg(feature = "arbitrary")]
fn is_subspace_basis<T: FiniteDimInnerSpace<Real = f64> + Display>(vs: &[T]) -> bool {
for i in 0 .. vs.len() {
// Basis elements must be normalized.
if !relative_eq!(vs[i].norm(), 1.0, epsilon = 1.0e-7) {
println!("Non-zero basis element norm: {}", vs[i].norm());
return false;
} }
fn $orthonormalization(vs: Vec<$Vector<f64>>) -> bool { for j in 0 .. i {
let mut basis = vs.clone(); // Basis elements must be orthogonal.
let subdim = $Vector::orthonormalize(&mut basis[..]); if !relative_eq!(vs[i].dot(&vs[j]), 0.0, epsilon = 1.0e-7) {
println!("Non-orthogonal basis elements: {} · {} = {}", vs[i], vs[j], vs[i].dot(&vs[j]));
if !is_subspace_basis(&basis[.. subdim]) { return false
return false;
} }
for mut e in vs {
for b in &basis[.. subdim] {
e -= e.dot(b) * b
}
// Any element of `e` must be a linear combination of the basis elements.
if !relative_eq!(e.norm(), 0.0, epsilon = 1.0e-7) {
println!("Orthonormalization; element decomposition failure: {}", e);
println!("... the non-zero norm is: {}", e.norm());
return false;
}
}
true
} }
} }
)*}
);
finite_dim_inner_space_test!( true
Vector1, orthonormal_subspace_basis1, orthonormalize1;
Vector2, orthonormal_subspace_basis2, orthonormalize2;
Vector3, orthonormal_subspace_basis3, orthonormalize3;
Vector4, orthonormal_subspace_basis4, orthonormalize4;
Vector5, orthonormal_subspace_basis5, orthonormalize5;
Vector6, orthonormal_subspace_basis6, orthonormalize6;
);
/*
*
* Helper functions.
*
*/
fn is_subspace_basis<T: FiniteDimInnerSpace<Real = f64> + Display>(vs: &[T]) -> bool {
for i in 0 .. vs.len() {
// Basis elements must be normalized.
if !relative_eq!(vs[i].norm(), 1.0, epsilon = 1.0e-7) {
println!("Non-zero basis element norm: {}", vs[i].norm());
return false;
}
for j in 0 .. i {
// Basis elements must be orthogonal.
if !relative_eq!(vs[i].dot(&vs[j]), 0.0, epsilon = 1.0e-7) {
println!("Non-orthogonal basis elements: {} · {} = {}", vs[i], vs[j], vs[i].dot(&vs[j]));
return false
}
}
} }
true
} }

View File

@ -1,4 +1,4 @@
use na::{Point3, Perspective3, Orthographic3}; use na::{Perspective3, Orthographic3};
#[test] #[test]
fn perspective_inverse() { fn perspective_inverse() {
@ -22,22 +22,26 @@ fn orthographic_inverse() {
#[cfg(feature = "arbitrary")] #[cfg(feature = "arbitrary")]
quickcheck!{ mod quickcheck_tests {
fn perspective_project_unproject(pt: Point3<f64>) -> bool { use na::{Point3, Perspective3, Orthographic3};
let proj = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0);
let projected = proj.project_point(&pt); quickcheck!{
let unprojected = proj.unproject_point(&projected); fn perspective_project_unproject(pt: Point3<f64>) -> bool {
let proj = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0);
relative_eq!(pt, unprojected, epsilon = 1.0e-7) let projected = proj.project_point(&pt);
} let unprojected = proj.unproject_point(&projected);
fn orthographic_project_unproject(pt: Point3<f64>) -> bool { relative_eq!(pt, unprojected, epsilon = 1.0e-7)
let proj = Orthographic3::new(1.0, 2.0, -3.0, -2.5, 10.0, 900.0); }
let projected = proj.project_point(&pt); fn orthographic_project_unproject(pt: Point3<f64>) -> bool {
let unprojected = proj.unproject_point(&projected); let proj = Orthographic3::new(1.0, 2.0, -3.0, -2.5, 10.0, 900.0);
relative_eq!(pt, unprojected, epsilon = 1.0e-7) let projected = proj.project_point(&pt);
let unprojected = proj.unproject_point(&projected);
relative_eq!(pt, unprojected, epsilon = 1.0e-7)
}
} }
} }

View File

@ -1,6 +1,4 @@
use std::f64; use na::{Vector2, Vector3};
use alga::general::Real;
use na::{self, Vector2, Vector3, Rotation2, Rotation3, Unit};
#[test] #[test]
fn angle_2() { fn angle_2() {
@ -19,187 +17,193 @@ fn angle_3() {
} }
#[cfg(feature = "arbitrary")] #[cfg(feature = "arbitrary")]
quickcheck!( mod quickcheck_tests {
/* use std::f64;
* use alga::general::Real;
* Euler angles. use na::{self, Vector2, Vector3, Rotation2, Rotation3, Unit};
*
*/
fn from_euler_angles(r: f64, p: f64, y: f64) -> bool {
let roll = Rotation3::from_euler_angles(r, 0.0, 0.0);
let pitch = Rotation3::from_euler_angles(0.0, p, 0.0);
let yaw = Rotation3::from_euler_angles(0.0, 0.0, y);
let rpy = Rotation3::from_euler_angles(r, p, y); quickcheck! {
/*
*
* Euler angles.
*
*/
fn from_euler_angles(r: f64, p: f64, y: f64) -> bool {
let roll = Rotation3::from_euler_angles(r, 0.0, 0.0);
let pitch = Rotation3::from_euler_angles(0.0, p, 0.0);
let yaw = Rotation3::from_euler_angles(0.0, 0.0, y);
roll[(0, 0)] == 1.0 && // rotation wrt. x axis. let rpy = Rotation3::from_euler_angles(r, p, y);
pitch[(1, 1)] == 1.0 && // rotation wrt. y axis.
yaw[(2, 2)] == 1.0 && // rotation wrt. z axis.
yaw * pitch * roll == rpy
}
fn to_euler_angles(r: f64, p: f64, y: f64) -> bool { roll[(0, 0)] == 1.0 && // rotation wrt. x axis.
let rpy = Rotation3::from_euler_angles(r, p, y); pitch[(1, 1)] == 1.0 && // rotation wrt. y axis.
let (roll, pitch, yaw) = rpy.to_euler_angles(); yaw[(2, 2)] == 1.0 && // rotation wrt. z axis.
relative_eq!(Rotation3::from_euler_angles(roll, pitch, yaw), rpy, epsilon = 1.0e-7) yaw * pitch * roll == rpy
}
fn to_euler_angles_gimble_lock(r: f64, y: f64) -> bool {
let pos = Rotation3::from_euler_angles(r, f64::frac_pi_2(), y);
let neg = Rotation3::from_euler_angles(r, -f64::frac_pi_2(), y);
let (pos_r, pos_p, pos_y) = pos.to_euler_angles();
let (neg_r, neg_p, neg_y) = neg.to_euler_angles();
relative_eq!(Rotation3::from_euler_angles(pos_r, pos_p, pos_y), pos, epsilon = 1.0e-7) &&
relative_eq!(Rotation3::from_euler_angles(neg_r, neg_p, neg_y), neg, epsilon = 1.0e-7)
}
/*
*
* Inversion is transposition.
*
*/
fn rotation_inv_3(a: Rotation3<f64>) -> bool {
let ta = a.transpose();
let ia = a.inverse();
ta == ia &&
relative_eq!(&ta * &a, Rotation3::identity(), epsilon = 1.0e-7) &&
relative_eq!(&ia * a, Rotation3::identity(), epsilon = 1.0e-7) &&
relative_eq!( a * &ta, Rotation3::identity(), epsilon = 1.0e-7) &&
relative_eq!( a * ia, Rotation3::identity(), epsilon = 1.0e-7)
}
fn rotation_inv_2(a: Rotation2<f64>) -> bool {
let ta = a.transpose();
let ia = a.inverse();
ta == ia &&
relative_eq!(&ta * &a, Rotation2::identity(), epsilon = 1.0e-7) &&
relative_eq!(&ia * a, Rotation2::identity(), epsilon = 1.0e-7) &&
relative_eq!( a * &ta, Rotation2::identity(), epsilon = 1.0e-7) &&
relative_eq!( a * ia, Rotation2::identity(), epsilon = 1.0e-7)
}
/*
*
* Angle between vectors.
*
*/
fn angle_is_commutative_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
a.angle(&b) == b.angle(&a)
}
fn angle_is_commutative_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
a.angle(&b) == b.angle(&a)
}
/*
*
* Rotation matrix between vectors.
*
*/
fn rotation_between_is_anticommutative_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
let rab = Rotation2::rotation_between(&a, &b);
let rba = Rotation2::rotation_between(&b, &a);
relative_eq!(rab * rba, Rotation2::identity())
}
fn rotation_between_is_anticommutative_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
let rots = (Rotation3::rotation_between(&a, &b), Rotation3::rotation_between(&b, &a));
if let (Some(rab), Some(rba)) = rots {
relative_eq!(rab * rba, Rotation3::identity(), epsilon = 1.0e-7)
} }
else {
true fn to_euler_angles(r: f64, p: f64, y: f64) -> bool {
let rpy = Rotation3::from_euler_angles(r, p, y);
let (roll, pitch, yaw) = rpy.to_euler_angles();
relative_eq!(Rotation3::from_euler_angles(roll, pitch, yaw), rpy, epsilon = 1.0e-7)
} }
}
fn rotation_between_is_identity(v2: Vector2<f64>, v3: Vector3<f64>) -> bool { fn to_euler_angles_gimble_lock(r: f64, y: f64) -> bool {
let vv2 = 3.42 * v2; let pos = Rotation3::from_euler_angles(r, f64::frac_pi_2(), y);
let vv3 = 4.23 * v3; let neg = Rotation3::from_euler_angles(r, -f64::frac_pi_2(), y);
let (pos_r, pos_p, pos_y) = pos.to_euler_angles();
relative_eq!(v2.angle(&vv2), 0.0, epsilon = 1.0e-7) && let (neg_r, neg_p, neg_y) = neg.to_euler_angles();
relative_eq!(v3.angle(&vv3), 0.0, epsilon = 1.0e-7) && relative_eq!(Rotation3::from_euler_angles(pos_r, pos_p, pos_y), pos, epsilon = 1.0e-7) &&
relative_eq!(Rotation2::rotation_between(&v2, &vv2), Rotation2::identity()) && relative_eq!(Rotation3::from_euler_angles(neg_r, neg_p, neg_y), neg, epsilon = 1.0e-7)
Rotation3::rotation_between(&v3, &vv3).unwrap() == Rotation3::identity()
}
fn rotation_between_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
if !relative_eq!(a.angle(&b), 0.0, epsilon = 1.0e-7) {
let r = Rotation2::rotation_between(&a, &b);
relative_eq!((r * a).angle(&b), 0.0, epsilon = 1.0e-7)
} }
else {
true /*
*
* Inversion is transposition.
*
*/
fn rotation_inv_3(a: Rotation3<f64>) -> bool {
let ta = a.transpose();
let ia = a.inverse();
ta == ia &&
relative_eq!(&ta * &a, Rotation3::identity(), epsilon = 1.0e-7) &&
relative_eq!(&ia * a, Rotation3::identity(), epsilon = 1.0e-7) &&
relative_eq!( a * &ta, Rotation3::identity(), epsilon = 1.0e-7) &&
relative_eq!( a * ia, Rotation3::identity(), epsilon = 1.0e-7)
} }
}
fn rotation_between_3(a: Vector3<f64>, b: Vector3<f64>) -> bool { fn rotation_inv_2(a: Rotation2<f64>) -> bool {
if !relative_eq!(a.angle(&b), 0.0, epsilon = 1.0e-7) { let ta = a.transpose();
let r = Rotation3::rotation_between(&a, &b).unwrap(); let ia = a.inverse();
relative_eq!((r * a).angle(&b), 0.0, epsilon = 1.0e-7)
ta == ia &&
relative_eq!(&ta * &a, Rotation2::identity(), epsilon = 1.0e-7) &&
relative_eq!(&ia * a, Rotation2::identity(), epsilon = 1.0e-7) &&
relative_eq!( a * &ta, Rotation2::identity(), epsilon = 1.0e-7) &&
relative_eq!( a * ia, Rotation2::identity(), epsilon = 1.0e-7)
} }
else {
true /*
*
* Angle between vectors.
*
*/
fn angle_is_commutative_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
a.angle(&b) == b.angle(&a)
}
fn angle_is_commutative_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
a.angle(&b) == b.angle(&a)
}
/*
*
* Rotation matrix between vectors.
*
*/
fn rotation_between_is_anticommutative_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
let rab = Rotation2::rotation_between(&a, &b);
let rba = Rotation2::rotation_between(&b, &a);
relative_eq!(rab * rba, Rotation2::identity())
}
fn rotation_between_is_anticommutative_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
let rots = (Rotation3::rotation_between(&a, &b), Rotation3::rotation_between(&b, &a));
if let (Some(rab), Some(rba)) = rots {
relative_eq!(rab * rba, Rotation3::identity(), epsilon = 1.0e-7)
}
else {
true
}
}
fn rotation_between_is_identity(v2: Vector2<f64>, v3: Vector3<f64>) -> bool {
let vv2 = 3.42 * v2;
let vv3 = 4.23 * v3;
relative_eq!(v2.angle(&vv2), 0.0, epsilon = 1.0e-7) &&
relative_eq!(v3.angle(&vv3), 0.0, epsilon = 1.0e-7) &&
relative_eq!(Rotation2::rotation_between(&v2, &vv2), Rotation2::identity()) &&
Rotation3::rotation_between(&v3, &vv3).unwrap() == Rotation3::identity()
}
fn rotation_between_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
if !relative_eq!(a.angle(&b), 0.0, epsilon = 1.0e-7) {
let r = Rotation2::rotation_between(&a, &b);
relative_eq!((r * a).angle(&b), 0.0, epsilon = 1.0e-7)
}
else {
true
}
}
fn rotation_between_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
if !relative_eq!(a.angle(&b), 0.0, epsilon = 1.0e-7) {
let r = Rotation3::rotation_between(&a, &b).unwrap();
relative_eq!((r * a).angle(&b), 0.0, epsilon = 1.0e-7)
}
else {
true
}
} }
}
/* /*
* *
* Rotation construction. * Rotation construction.
* *
*/ */
fn new_rotation_2(angle: f64) -> bool { fn new_rotation_2(angle: f64) -> bool {
let r = Rotation2::new(angle); let r = Rotation2::new(angle);
let angle = na::wrap(angle, -f64::pi(), f64::pi());
relative_eq!(r.angle(), angle, epsilon = 1.0e-7)
}
fn new_rotation_3(axisangle: Vector3<f64>) -> bool {
let r = Rotation3::new(axisangle);
if let Some((axis, angle)) = Unit::try_new_and_get(axisangle, 0.0) {
let angle = na::wrap(angle, -f64::pi(), f64::pi()); let angle = na::wrap(angle, -f64::pi(), f64::pi());
(relative_eq!(r.angle(), angle, epsilon = 1.0e-7) && relative_eq!(r.angle(), angle, epsilon = 1.0e-7)
relative_eq!(r.axis().unwrap(), axis, epsilon = 1.0e-7)) ||
(relative_eq!(r.angle(), -angle, epsilon = 1.0e-7) &&
relative_eq!(r.axis().unwrap(), -axis, epsilon = 1.0e-7))
} }
else {
r == Rotation3::identity() fn new_rotation_3(axisangle: Vector3<f64>) -> bool {
let r = Rotation3::new(axisangle);
if let Some((axis, angle)) = Unit::try_new_and_get(axisangle, 0.0) {
let angle = na::wrap(angle, -f64::pi(), f64::pi());
(relative_eq!(r.angle(), angle, epsilon = 1.0e-7) &&
relative_eq!(r.axis().unwrap(), axis, epsilon = 1.0e-7)) ||
(relative_eq!(r.angle(), -angle, epsilon = 1.0e-7) &&
relative_eq!(r.axis().unwrap(), -axis, epsilon = 1.0e-7))
}
else {
r == Rotation3::identity()
}
} }
}
/* /*
* *
* Rotation pow. * Rotation pow.
* *
*/ */
fn powf_rotation_2(angle: f64, pow: f64) -> bool { fn powf_rotation_2(angle: f64, pow: f64) -> bool {
let r = Rotation2::new(angle).powf(pow); let r = Rotation2::new(angle).powf(pow);
let angle = na::wrap(angle, -f64::pi(), f64::pi()); let angle = na::wrap(angle, -f64::pi(), f64::pi());
let pangle = na::wrap(angle * pow, -f64::pi(), f64::pi());
relative_eq!(r.angle(), pangle, epsilon = 1.0e-7)
}
fn powf_rotation_3(axisangle: Vector3<f64>, pow: f64) -> bool {
let r = Rotation3::new(axisangle).powf(pow);
if let Some((axis, angle)) = Unit::try_new_and_get(axisangle, 0.0) {
let angle = na::wrap(angle, -f64::pi(), f64::pi());
let pangle = na::wrap(angle * pow, -f64::pi(), f64::pi()); let pangle = na::wrap(angle * pow, -f64::pi(), f64::pi());
relative_eq!(r.angle(), pangle, epsilon = 1.0e-7)
(relative_eq!(r.angle(), pangle, epsilon = 1.0e-7) &&
relative_eq!(r.axis().unwrap(), axis, epsilon = 1.0e-7)) ||
(relative_eq!(r.angle(), -pangle, epsilon = 1.0e-7) &&
relative_eq!(r.axis().unwrap(), -axis, epsilon = 1.0e-7))
} }
else {
r == Rotation3::identity() fn powf_rotation_3(axisangle: Vector3<f64>, pow: f64) -> bool {
let r = Rotation3::new(axisangle).powf(pow);
if let Some((axis, angle)) = Unit::try_new_and_get(axisangle, 0.0) {
let angle = na::wrap(angle, -f64::pi(), f64::pi());
let pangle = na::wrap(angle * pow, -f64::pi(), f64::pi());
(relative_eq!(r.angle(), pangle, epsilon = 1.0e-7) &&
relative_eq!(r.axis().unwrap(), axis, epsilon = 1.0e-7)) ||
(relative_eq!(r.angle(), -pangle, epsilon = 1.0e-7) &&
relative_eq!(r.axis().unwrap(), -axis, epsilon = 1.0e-7))
}
else {
r == Rotation3::identity()
}
} }
} }
); }

View File

@ -1,9 +1,9 @@
#![cfg(feature = "arbitrary")]
#![allow(non_snake_case)] #![allow(non_snake_case)]
use alga::linear::{Transformation, ProjectiveTransformation}; use alga::linear::{Transformation, ProjectiveTransformation};
use na::{Vector3, Point3, Similarity3, Translation3, Isometry3, UnitQuaternion}; use na::{Vector3, Point3, Similarity3, Translation3, Isometry3, UnitQuaternion};
#[cfg(feature = "arbitrary")]
quickcheck!( quickcheck!(
fn inverse_is_identity(i: Similarity3<f64>, p: Point3<f64>, v: Vector3<f64>) -> bool { fn inverse_is_identity(i: Similarity3<f64>, p: Point3<f64>, v: Vector3<f64>) -> bool {
let ii = i.inverse(); let ii = i.inverse();

View File

@ -1,9 +1,8 @@
#![cfg(feature = "arbitrary")]
#![allow(non_snake_case)] #![allow(non_snake_case)]
use na::{Unit, UnitComplex, Vector2, Point2, Rotation2}; use na::{Unit, UnitComplex, Vector2, Point2, Rotation2};
#[cfg(feature = "arbitrary")]
quickcheck!( quickcheck!(
/* /*

View File

@ -1,9 +1,10 @@
#![cfg(feature = "arbitrary")]
use std::cmp; use std::cmp;
use na::{DMatrix, Matrix4}; use na::{DMatrix, Matrix4};
use na::balancing; use na::balancing;
#[cfg(feature = "arbitrary")]
quickcheck! { quickcheck! {
fn balancing_parlett_reinsch(n: usize) -> bool { fn balancing_parlett_reinsch(n: usize) -> bool {
let n = cmp::min(n, 10); let n = cmp::min(n, 10);

View File

@ -1,7 +1,7 @@
#![cfg(feature = "arbitrary")]
use na::{DMatrix, Matrix2, Matrix4, Matrix5x3, Matrix3x5}; use na::{DMatrix, Matrix2, Matrix4, Matrix5x3, Matrix3x5};
#[cfg(feature = "arbitrary")]
quickcheck! { quickcheck! {
fn bidiagonal(m: DMatrix<f64>) -> bool { fn bidiagonal(m: DMatrix<f64>) -> bool {
if m.len() == 0 { if m.len() == 0 {

View File

@ -1,58 +1,61 @@
use std::cmp; use na::DMatrix;
use na::{DMatrix, Matrix2, Matrix3, Matrix4};
#[cfg(feature = "arbitrary")] #[cfg(feature = "arbitrary")]
quickcheck! { mod quickcheck_tests {
fn symmetric_eigen(n: usize) -> bool { use std::cmp;
let n = cmp::max(1, cmp::min(n, 10)); use na::{DMatrix, Matrix2, Matrix3, Matrix4};
let m = DMatrix::<f64>::new_random(n, n);
let eig = m.clone().symmetric_eigen();
let recomp = eig.recompose();
println!("{}{}", m.lower_triangle(), recomp.lower_triangle()); 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 = m.clone().symmetric_eigen();
let recomp = eig.recompose();
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5) println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
}
fn symmetric_eigen_singular(n: usize) -> bool { relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
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 = m.clone().symmetric_eigen();
let recomp = eig.recompose();
println!("{}{}", m.lower_triangle(), recomp.lower_triangle()); 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 = m.clone().symmetric_eigen();
let recomp = eig.recompose();
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5) println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
}
fn symmetric_eigen_static_square_4x4(m: Matrix4<f64>) -> bool { relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
let eig = m.symmetric_eigen(); }
let recomp = eig.recompose();
println!("{}{}", m.lower_triangle(), recomp.lower_triangle()); fn symmetric_eigen_static_square_4x4(m: Matrix4<f64>) -> bool {
let eig = m.symmetric_eigen();
let recomp = eig.recompose();
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5) println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
}
fn symmetric_eigen_static_square_3x3(m: Matrix3<f64>) -> bool { relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
let eig = m.symmetric_eigen(); }
let recomp = eig.recompose();
println!("{}{}", m.lower_triangle(), recomp.lower_triangle()); fn symmetric_eigen_static_square_3x3(m: Matrix3<f64>) -> bool {
let eig = m.symmetric_eigen();
let recomp = eig.recompose();
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5) println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
}
fn symmetric_eigen_static_square_2x2(m: Matrix2<f64>) -> bool { relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
let eig = m.symmetric_eigen(); }
let recomp = eig.recompose();
println!("{}{}", m.lower_triangle(), recomp.lower_triangle()); fn symmetric_eigen_static_square_2x2(m: Matrix2<f64>) -> bool {
let eig = m.symmetric_eigen();
let recomp = eig.recompose();
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5) println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
}
} }
} }

View File

@ -1,7 +1,4 @@
use std::cmp; use na::Matrix3;
use na::{DMatrix, Matrix3, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5,
DVector, Vector4};
#[test] #[test]
fn full_piv_lu_simple() { fn full_piv_lu_simple() {
@ -42,114 +39,119 @@ fn full_piv_lu_simple_with_pivot() {
} }
#[cfg(feature = "arbitrary")] #[cfg(feature = "arbitrary")]
quickcheck! { mod quickcheck_tests {
fn full_piv_lu(m: DMatrix<f64>) -> bool { use std::cmp;
let mut m = m; use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, DVector, Vector4};
if m.len() == 0 {
m = DMatrix::new_random(1, 1);
}
let lu = m.clone().full_piv_lu(); quickcheck! {
let (p, l, u, q) = lu.unpack(); fn full_piv_lu(m: DMatrix<f64>) -> bool {
let mut lu = l * u; let mut m = m;
p.inv_permute_rows(&mut lu); if m.len() == 0 {
q.inv_permute_columns(&mut lu); m = DMatrix::new_random(1, 1);
}
relative_eq!(m, lu, epsilon = 1.0e-7)
}
fn full_piv_lu_static_3_5(m: Matrix3x5<f64>) -> bool {
let lu = m.full_piv_lu();
let (p, l, u, q) = lu.unpack();
let mut lu = l * u;
p.inv_permute_rows(&mut lu);
q.inv_permute_columns(&mut lu);
relative_eq!(m, lu, epsilon = 1.0e-7)
}
fn full_piv_lu_static_5_3(m: Matrix5x3<f64>) -> bool {
let lu = m.full_piv_lu();
let (p, l, u, q) = lu.unpack();
let mut lu = l * u;
p.inv_permute_rows(&mut lu);
q.inv_permute_columns(&mut lu);
relative_eq!(m, lu, epsilon = 1.0e-7)
}
fn full_piv_lu_static_square(m: Matrix4<f64>) -> bool {
let lu = m.full_piv_lu();
let (p, l, u, q) = lu.unpack();
let mut lu = l * u;
p.inv_permute_rows(&mut lu);
q.inv_permute_columns(&mut lu);
relative_eq!(m, lu, epsilon = 1.0e-7)
}
fn full_piv_lu_solve(n: usize, nb: usize) -> bool {
if n != 0 && nb != 0 {
let n = cmp::min(n, 50); // To avoid slowing down the test too much.
let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
let m = DMatrix::<f64>::new_random(n, n);
let lu = m.clone().full_piv_lu(); let lu = m.clone().full_piv_lu();
let b1 = DVector::new_random(n); let (p, l, u, q) = lu.unpack();
let b2 = DMatrix::new_random(n, nb); let mut lu = l * u;
p.inv_permute_rows(&mut lu);
q.inv_permute_columns(&mut lu);
let sol1 = lu.solve(&b1); relative_eq!(m, lu, epsilon = 1.0e-7)
let sol2 = lu.solve(&b2);
return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
} }
return true; fn full_piv_lu_static_3_5(m: Matrix3x5<f64>) -> bool {
} let lu = m.full_piv_lu();
let (p, l, u, q) = lu.unpack();
let mut lu = l * u;
p.inv_permute_rows(&mut lu);
q.inv_permute_columns(&mut lu);
fn full_piv_lu_solve_static(m: Matrix4<f64>) -> bool { relative_eq!(m, lu, epsilon = 1.0e-7)
let lu = m.full_piv_lu(); }
let b1 = Vector4::new_random();
let b2 = Matrix4x3::new_random();
let sol1 = lu.solve(&b1); fn full_piv_lu_static_5_3(m: Matrix5x3<f64>) -> bool {
let sol2 = lu.solve(&b2); let lu = m.full_piv_lu();
let (p, l, u, q) = lu.unpack();
let mut lu = l * u;
p.inv_permute_rows(&mut lu);
q.inv_permute_columns(&mut lu);
return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) && relative_eq!(m, lu, epsilon = 1.0e-7)
(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6)) }
}
fn full_piv_lu_inverse(n: usize) -> bool { fn full_piv_lu_static_square(m: Matrix4<f64>) -> bool {
let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much. let lu = m.full_piv_lu();
let m = DMatrix::<f64>::new_random(n, n); let (p, l, u, q) = lu.unpack();
let mut lu = l * u;
p.inv_permute_rows(&mut lu);
q.inv_permute_columns(&mut lu);
let mut l = m.lower_triangle(); relative_eq!(m, lu, epsilon = 1.0e-7)
let mut u = m.upper_triangle(); }
// Ensure the matrix is well conditioned for inversion. fn full_piv_lu_solve(n: usize, nb: usize) -> bool {
l.fill_diagonal(1.0); if n != 0 && nb != 0 {
u.fill_diagonal(1.0); let n = cmp::min(n, 50); // To avoid slowing down the test too much.
let m = l * u; let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
let m = DMatrix::<f64>::new_random(n, n);
let m1 = m.clone().full_piv_lu().try_inverse().unwrap(); let lu = m.clone().full_piv_lu();
let id1 = &m * &m1; let b1 = DVector::new_random(n);
let id2 = &m1 * &m; let b2 = DMatrix::new_random(n, nb);
return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5); let sol1 = lu.solve(&b1);
} let sol2 = lu.solve(&b2);
fn full_piv_lu_inverse_static(m: Matrix4<f64>) -> bool { return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
let lu = m.full_piv_lu(); (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
}
if let Some(m1) = lu.try_inverse() { return true;
}
fn full_piv_lu_solve_static(m: Matrix4<f64>) -> bool {
let lu = m.full_piv_lu();
let b1 = Vector4::new_random();
let b2 = Matrix4x3::new_random();
let sol1 = lu.solve(&b1);
let sol2 = lu.solve(&b2);
return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
}
fn full_piv_lu_inverse(n: usize) -> bool {
let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
let m = DMatrix::<f64>::new_random(n, n);
let mut l = m.lower_triangle();
let mut u = m.upper_triangle();
// Ensure the matrix is well conditioned for inversion.
l.fill_diagonal(1.0);
u.fill_diagonal(1.0);
let m = l * u;
let m1 = m.clone().full_piv_lu().try_inverse().unwrap();
let id1 = &m * &m1; let id1 = &m * &m1;
let id2 = &m1 * &m; let id2 = &m1 * &m;
id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5) return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5);
} }
else {
true fn full_piv_lu_inverse_static(m: Matrix4<f64>) -> bool {
let lu = m.full_piv_lu();
if let Some(m1) = lu.try_inverse() {
let id1 = &m * &m1;
let id2 = &m1 * &m;
id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
}
else {
true
}
} }
} }
} }

View File

@ -1,7 +1,8 @@
#![cfg(feature = "arbitrary")]
use std::cmp; use std::cmp;
use na::{DMatrix, Matrix2, Matrix4}; use na::{DMatrix, Matrix2, Matrix4};
#[test] #[test]
fn hessenberg_simple() { fn hessenberg_simple() {
let m = Matrix2::new(1.0, 0.0, let m = Matrix2::new(1.0, 0.0,
@ -11,8 +12,6 @@ fn hessenberg_simple() {
assert!(relative_eq!(m, p * h * p.transpose(), epsilon = 1.0e-7)) assert!(relative_eq!(m, p * h * p.transpose(), epsilon = 1.0e-7))
} }
#[cfg(feature = "arbitrary")]
quickcheck! { quickcheck! {
fn hessenberg(n: usize) -> bool { fn hessenberg(n: usize) -> bool {
let n = cmp::max(1, cmp::min(n, 50)); let n = cmp::max(1, cmp::min(n, 50));

View File

@ -1,7 +1,4 @@
use std::cmp; use na::Matrix3;
use na::{DMatrix, Matrix3, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5,
DVector, Vector4};
#[test] #[test]
fn lu_simple() { fn lu_simple() {
@ -40,110 +37,115 @@ fn lu_simple_with_pivot() {
} }
#[cfg(feature = "arbitrary")] #[cfg(feature = "arbitrary")]
quickcheck! { mod quickcheck_tests {
fn lu(m: DMatrix<f64>) -> bool { use std::cmp;
let mut m = m; use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, DVector, Vector4};
if m.len() == 0 {
m = DMatrix::new_random(1, 1);
}
let lu = m.clone().lu(); quickcheck! {
let (p, l, u) = lu.unpack(); fn lu(m: DMatrix<f64>) -> bool {
let mut lu = l * u; let mut m = m;
p.inv_permute_rows(&mut lu); if m.len() == 0 {
m = DMatrix::new_random(1, 1);
relative_eq!(m, lu, epsilon = 1.0e-7) }
}
fn lu_static_3_5(m: Matrix3x5<f64>) -> bool {
let lu = m.lu();
let (p, l, u) = lu.unpack();
let mut lu = l * u;
p.inv_permute_rows(&mut lu);
relative_eq!(m, lu, epsilon = 1.0e-7)
}
fn lu_static_5_3(m: Matrix5x3<f64>) -> bool {
let lu = m.lu();
let (p, l, u) = lu.unpack();
let mut lu = l * u;
p.inv_permute_rows(&mut lu);
relative_eq!(m, lu, epsilon = 1.0e-7)
}
fn lu_static_square(m: Matrix4<f64>) -> bool {
let lu = m.lu();
let (p, l, u) = lu.unpack();
let mut lu = l * u;
p.inv_permute_rows(&mut lu);
relative_eq!(m, lu, epsilon = 1.0e-7)
}
fn lu_solve(n: usize, nb: usize) -> bool {
if n != 0 && nb != 0 {
let n = cmp::min(n, 50); // To avoid slowing down the test too much.
let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
let m = DMatrix::<f64>::new_random(n, n);
let lu = m.clone().lu(); let lu = m.clone().lu();
let b1 = DVector::new_random(n); let (p, l, u) = lu.unpack();
let b2 = DMatrix::new_random(n, nb); let mut lu = l * u;
p.inv_permute_rows(&mut lu);
let sol1 = lu.solve(&b1); relative_eq!(m, lu, epsilon = 1.0e-7)
let sol2 = lu.solve(&b2);
return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
} }
return true; fn lu_static_3_5(m: Matrix3x5<f64>) -> bool {
} let lu = m.lu();
let (p, l, u) = lu.unpack();
let mut lu = l * u;
p.inv_permute_rows(&mut lu);
fn lu_solve_static(m: Matrix4<f64>) -> bool { relative_eq!(m, lu, epsilon = 1.0e-7)
let lu = m.lu(); }
let b1 = Vector4::new_random();
let b2 = Matrix4x3::new_random();
let sol1 = lu.solve(&b1); fn lu_static_5_3(m: Matrix5x3<f64>) -> bool {
let sol2 = lu.solve(&b2); let lu = m.lu();
let (p, l, u) = lu.unpack();
let mut lu = l * u;
p.inv_permute_rows(&mut lu);
return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) && relative_eq!(m, lu, epsilon = 1.0e-7)
(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6)) }
}
fn lu_inverse(n: usize) -> bool { fn lu_static_square(m: Matrix4<f64>) -> bool {
let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much. let lu = m.lu();
let m = DMatrix::<f64>::new_random(n, n); let (p, l, u) = lu.unpack();
let mut lu = l * u;
p.inv_permute_rows(&mut lu);
let mut l = m.lower_triangle(); relative_eq!(m, lu, epsilon = 1.0e-7)
let mut u = m.upper_triangle(); }
// Ensure the matrix is well conditioned for inversion. fn lu_solve(n: usize, nb: usize) -> bool {
l.fill_diagonal(1.0); if n != 0 && nb != 0 {
u.fill_diagonal(1.0); let n = cmp::min(n, 50); // To avoid slowing down the test too much.
let m = l * u; let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
let m = DMatrix::<f64>::new_random(n, n);
let m1 = m.clone().lu().try_inverse().unwrap(); let lu = m.clone().lu();
let id1 = &m * &m1; let b1 = DVector::new_random(n);
let id2 = &m1 * &m; let b2 = DMatrix::new_random(n, nb);
return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5); let sol1 = lu.solve(&b1);
} let sol2 = lu.solve(&b2);
fn lu_inverse_static(m: Matrix4<f64>) -> bool { return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
let lu = m.lu(); (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
}
if let Some(m1) = lu.try_inverse() { return true;
}
fn lu_solve_static(m: Matrix4<f64>) -> bool {
let lu = m.lu();
let b1 = Vector4::new_random();
let b2 = Matrix4x3::new_random();
let sol1 = lu.solve(&b1);
let sol2 = lu.solve(&b2);
return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
}
fn lu_inverse(n: usize) -> bool {
let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
let m = DMatrix::<f64>::new_random(n, n);
let mut l = m.lower_triangle();
let mut u = m.upper_triangle();
// Ensure the matrix is well conditioned for inversion.
l.fill_diagonal(1.0);
u.fill_diagonal(1.0);
let m = l * u;
let m1 = m.clone().lu().try_inverse().unwrap();
let id1 = &m * &m1; let id1 = &m * &m1;
let id2 = &m1 * &m; let id2 = &m1 * &m;
id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5) return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5);
} }
else {
true fn lu_inverse_static(m: Matrix4<f64>) -> bool {
let lu = m.lu();
if let Some(m1) = lu.try_inverse() {
let id1 = &m * &m1;
let id2 = &m1 * &m;
id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
}
else {
true
}
} }
} }
} }

View File

@ -1,8 +1,9 @@
#![cfg(feature = "arbitrary")]
use std::cmp; use std::cmp;
use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5,
DVector, Vector4}; DVector, Vector4};
#[cfg(feature = "arbitrary")]
quickcheck! { quickcheck! {
fn qr(m: DMatrix<f64>) -> bool { fn qr(m: DMatrix<f64>) -> bool {
let qr = m.clone().qr(); let qr = m.clone().qr();

View File

@ -1,5 +1,4 @@
use std::cmp; use na::{DMatrix, Matrix3, Matrix4};
use na::{DMatrix, Matrix2, Matrix3, Matrix4};
#[test] #[test]
@ -15,49 +14,54 @@ fn schur_simpl_mat3() {
} }
#[cfg(feature = "arbitrary")] #[cfg(feature = "arbitrary")]
quickcheck! { mod quickcheck_tests {
fn schur(n: usize) -> bool { use std::cmp;
let n = cmp::max(1, cmp::min(n, 10)); use na::{DMatrix, Matrix2, Matrix3, Matrix4};
let m = DMatrix::<f64>::new_random(n, n);
let (vecs, vals) = m.clone().real_schur().unpack(); quickcheck! {
fn schur(n: usize) -> bool {
let n = cmp::max(1, cmp::min(n, 10));
let m = DMatrix::<f64>::new_random(n, n);
if !relative_eq!(&vecs * &vals * vecs.transpose(), m, epsilon = 1.0e-7) { let (vecs, vals) = m.clone().real_schur().unpack();
println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
if !relative_eq!(&vecs * &vals * vecs.transpose(), m, epsilon = 1.0e-7) {
println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
}
relative_eq!(&vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7)
} }
relative_eq!(&vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7) fn schur_static_mat2(m: Matrix2<f64>) -> bool {
} let (vecs, vals) = m.clone().real_schur().unpack();
fn schur_static_mat2(m: Matrix2<f64>) -> bool { let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7);
let (vecs, vals) = m.clone().real_schur().unpack(); if !ok {
println!("{:.5}{:.5}", vecs, vals);
let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7); println!("Reconstruction:{}{}", m, &vecs * &vals * vecs.transpose());
if !ok { }
println!("{:.5}{:.5}", vecs, vals); ok
println!("Reconstruction:{}{}", m, &vecs * &vals * vecs.transpose());
} }
ok
}
fn schur_static_mat3(m: Matrix3<f64>) -> bool { fn schur_static_mat3(m: Matrix3<f64>) -> bool {
let (vecs, vals) = m.clone().real_schur().unpack(); let (vecs, vals) = m.clone().real_schur().unpack();
let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7); let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7);
if !ok { if !ok {
println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose()); println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
}
ok
} }
ok
}
fn schur_static_mat4(m: Matrix4<f64>) -> bool { fn schur_static_mat4(m: Matrix4<f64>) -> bool {
let (vecs, vals) = m.clone().real_schur().unpack(); let (vecs, vals) = m.clone().real_schur().unpack();
let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7); let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7);
if !ok { if !ok {
println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose()); println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
}
ok
} }
ok
} }
} }

View File

@ -1,3 +1,5 @@
#![cfg(feature = "arbitrary")]
use na::{Matrix4, Matrix4x5}; use na::{Matrix4, Matrix4x5};
fn unzero_diagonal(a: &mut Matrix4<f64>) { fn unzero_diagonal(a: &mut Matrix4<f64>) {
@ -8,7 +10,6 @@ fn unzero_diagonal(a: &mut Matrix4<f64>) {
} }
} }
#[cfg(feature = "arbitrary")]
quickcheck! { quickcheck! {
fn solve_lower_triangular(a: Matrix4<f64>, b: Matrix4x5<f64>) -> bool { fn solve_lower_triangular(a: Matrix4<f64>, b: Matrix4x5<f64>) -> bool {
let mut a = a; let mut a = a;

View File

@ -1,142 +1,143 @@
use std::cmp; use na::{DMatrix, Matrix6};
use na::{DMatrix, Matrix2, Matrix3, Matrix4, Matrix6, Matrix5x2, Matrix5x3, Matrix2x5, Matrix3x5,
DVector};
#[cfg(feature = "arbitrary")] #[cfg(feature = "arbitrary")]
quickcheck! { mod quickcheck_tests {
fn svd(m: DMatrix<f64>) -> bool { use std::cmp;
if m.len() > 0 { use na::{DMatrix, Matrix2, Matrix3, Matrix4, Matrix5x2, Matrix5x3, Matrix2x5, Matrix3x5, DVector};
let svd = m.clone().svd(true, true);
let recomp_m = svd.clone().recompose();
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = DMatrix::from_diagonal(&s);
println!("{}{}", &m, &u * &ds * &v_t); quickcheck! {
fn svd(m: DMatrix<f64>) -> bool {
if m.len() > 0 {
let svd = m.clone().svd(true, true);
let recomp_m = svd.clone().recompose();
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = DMatrix::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) && println!("{}{}", &m, &u * &ds * &v_t);
relative_eq!(&u * ds * &v_t, recomp_m, epsilon = 1.0e-5) &&
relative_eq!(m, recomp_m, epsilon = 1.0e-5)
}
else {
true
}
}
fn svd_static_5_3(m: Matrix5x3<f64>) -> bool { s.iter().all(|e| *e >= 0.0) &&
let svd = m.svd(true, true); relative_eq!(&u * ds * &v_t, recomp_m, epsilon = 1.0e-5) &&
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap()); relative_eq!(m, recomp_m, epsilon = 1.0e-5)
let ds = Matrix3::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5)
}
fn svd_static_5_2(m: Matrix5x2<f64>) -> bool {
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = Matrix2::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5)
}
fn svd_static_3_5(m: Matrix3x5<f64>) -> bool {
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = Matrix3::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
}
fn svd_static_2_5(m: Matrix2x5<f64>) -> bool {
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = Matrix2::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
}
fn svd_static_square(m: Matrix4<f64>) -> bool {
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = Matrix4::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5)
}
fn svd_static_square_2x2(m: Matrix2<f64>) -> bool {
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = Matrix2::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5)
}
fn svd_pseudo_inverse(m: DMatrix<f64>) -> bool {
if m.len() > 0 {
let svd = m.clone().svd(true, true);
let pinv = svd.pseudo_inverse(1.0e-10);
if m.nrows() > m.ncols() {
println!("{}", &pinv * &m);
(pinv * m).is_identity(1.0e-5)
} }
else { else {
println!("{}", &m * &pinv); true
(m * pinv).is_identity(1.0e-5)
} }
} }
else {
fn svd_static_5_3(m: Matrix5x3<f64>) -> bool {
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = Matrix3::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5)
}
fn svd_static_5_2(m: Matrix5x2<f64>) -> bool {
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = Matrix2::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5)
}
fn svd_static_3_5(m: Matrix3x5<f64>) -> bool {
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = Matrix3::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
}
fn svd_static_2_5(m: Matrix2x5<f64>) -> bool {
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = Matrix2::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
}
fn svd_static_square(m: Matrix4<f64>) -> bool {
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = Matrix4::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5)
}
fn svd_static_square_2x2(m: Matrix2<f64>) -> bool {
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = Matrix2::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5)
}
fn svd_pseudo_inverse(m: DMatrix<f64>) -> bool {
if m.len() > 0 {
let svd = m.clone().svd(true, true);
let pinv = svd.pseudo_inverse(1.0e-10);
if m.nrows() > m.ncols() {
println!("{}", &pinv * &m);
(pinv * m).is_identity(1.0e-5)
}
else {
println!("{}", &m * &pinv);
(m * pinv).is_identity(1.0e-5)
}
}
else {
true
}
}
fn svd_solve(n: usize, nb: usize) -> bool {
let n = cmp::max(1, cmp::min(n, 10));
let nb = cmp::min(nb, 10);
let m = DMatrix::<f64>::new_random(n, n);
let svd = m.clone().svd(true, true);
if svd.rank(1.0e-7) == n {
let b1 = DVector::new_random(n);
let b2 = DMatrix::new_random(n, nb);
let sol1 = svd.solve(&b1, 1.0e-7);
let sol2 = svd.solve(&b2, 1.0e-7);
let recomp = svd.recompose();
if !relative_eq!(m, recomp, epsilon = 1.0e-6) {
println!("{}{}", m, recomp);
}
if !relative_eq!(&m * &sol1, b1, epsilon = 1.0e-6) {
println!("Problem 1: {:.6}{:.6}", b1, &m * sol1);
return false;
}
if !relative_eq!(&m * &sol2, b2, epsilon = 1.0e-6) {
println!("Problem 2: {:.6}{:.6}", b2, &m * sol2);
return false;
}
}
true true
} }
} }
fn svd_solve(n: usize, nb: usize) -> bool {
let n = cmp::max(1, cmp::min(n, 10));
let nb = cmp::min(nb, 10);
let m = DMatrix::<f64>::new_random(n, n);
let svd = m.clone().svd(true, true);
if svd.rank(1.0e-7) == n {
let b1 = DVector::new_random(n);
let b2 = DMatrix::new_random(n, nb);
let sol1 = svd.solve(&b1, 1.0e-7);
let sol2 = svd.solve(&b2, 1.0e-7);
let recomp = svd.recompose();
if !relative_eq!(m, recomp, epsilon = 1.0e-6) {
println!("{}{}", m, recomp);
}
if !relative_eq!(&m * &sol1, b1, epsilon = 1.0e-6) {
println!("Problem 1: {:.6}{:.6}", b1, &m * sol1);
return false;
}
if !relative_eq!(&m * &sol2, b2, epsilon = 1.0e-6) {
println!("Problem 2: {:.6}{:.6}", b2, &m * sol2);
return false;
}
}
true
}
} }
// Test proposed on the issue #176 of rulinalg. // Test proposed on the issue #176 of rulinalg.

View File

@ -1,9 +1,9 @@
#![cfg(feature = "arbitrary")]
use std::cmp; use std::cmp;
use na::{DMatrix, Matrix2, Matrix4}; use na::{DMatrix, Matrix2, Matrix4};
#[cfg(feature = "arbitrary")]
quickcheck! { quickcheck! {
fn symm_tridiagonal(n: usize) -> bool { fn symm_tridiagonal(n: usize) -> bool {
let n = cmp::max(1, cmp::min(n, 50)); let n = cmp::max(1, cmp::min(n, 50));