From bba1993e58f9d03a16b6161dc1aab2381485b9bf Mon Sep 17 00:00:00 2001 From: Eduard Bopp Date: Wed, 17 Jan 2018 16:48:47 +0100 Subject: [PATCH] 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. --- tests/core/blas.rs | 5 +- tests/core/matrix.rs | 158 ++++++++-------- tests/geometry/projection.rs | 30 +-- tests/geometry/rotation.rs | 332 +++++++++++++++++---------------- tests/geometry/similarity.rs | 2 +- tests/geometry/unit_complex.rs | 3 +- tests/linalg/balancing.rs | 3 +- tests/linalg/bidiagonal.rs | 4 +- tests/linalg/eigen.rs | 81 ++++---- tests/linalg/full_piv_lu.rs | 186 +++++++++--------- tests/linalg/hessenberg.rs | 5 +- tests/linalg/lu.rs | 178 +++++++++--------- tests/linalg/qr.rs | 3 +- tests/linalg/real_schur.rs | 70 +++---- tests/linalg/solve.rs | 3 +- tests/linalg/svd.rs | 255 ++++++++++++------------- tests/linalg/tridiagonal.rs | 4 +- 17 files changed, 673 insertions(+), 649 deletions(-) diff --git a/tests/core/blas.rs b/tests/core/blas.rs index 57fab2c6..82239af0 100644 --- a/tests/core/blas.rs +++ b/tests/core/blas.rs @@ -1,9 +1,8 @@ -use std::cmp; +#![cfg(feature = "arbitrary")] +use std::cmp; use na::{DVector, DMatrix}; - -#[cfg(feature = "arbitrary")] quickcheck! { /* * diff --git a/tests/core/matrix.rs b/tests/core/matrix.rs index 2ebab858..996d1d64 100644 --- a/tests/core/matrix.rs +++ b/tests/core/matrix.rs @@ -1,8 +1,5 @@ use num::{Zero, One}; use num::Float; -use std::fmt::Display; - -use alga::linear::FiniteDimInnerSpace; use na::{self, DVector, DMatrix, @@ -806,93 +803,100 @@ mod normalization_tests { } } +#[cfg(feature = "arbitrary")] // FIXME: move this to alga ? -macro_rules! finite_dim_inner_space_test( - ($($Vector: ident, $orthonormal_subspace: ident, $orthonormalization: ident);* $(;)*) => {$( - #[cfg(feature = "arbitrary")] - quickcheck!{ - fn $orthonormal_subspace(vs: Vec<$Vector>) -> 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 } - ); +mod finite_dim_inner_space_tests { + use super::*; + use std::fmt::Display; + use alga::linear::FiniteDimInnerSpace; - if !is_subspace_basis(&ortho_basis[..]) { - return false; + macro_rules! finite_dim_inner_space_test( + ($($Vector: ident, $orthonormal_subspace: ident, $orthonormalization: ident);* $(;)*) => {$( + quickcheck!{ + fn $orthonormal_subspace(vs: Vec<$Vector>) -> 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 { - 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)); + fn $orthonormalization(vs: Vec<$Vector>) -> bool { + let mut basis = vs.clone(); + let subdim = $Vector::orthonormalize(&mut basis[..]); + + 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; } } - } - 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 + 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>) -> bool { - let mut basis = vs.clone(); - let subdim = $Vector::orthonormalize(&mut basis[..]); - - if !is_subspace_basis(&basis[.. subdim]) { - 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 } - - 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!( - 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 + 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 } - - true } diff --git a/tests/geometry/projection.rs b/tests/geometry/projection.rs index 84f3cf70..7be0d49d 100644 --- a/tests/geometry/projection.rs +++ b/tests/geometry/projection.rs @@ -1,4 +1,4 @@ -use na::{Point3, Perspective3, Orthographic3}; +use na::{Perspective3, Orthographic3}; #[test] fn perspective_inverse() { @@ -22,22 +22,26 @@ fn orthographic_inverse() { #[cfg(feature = "arbitrary")] -quickcheck!{ - fn perspective_project_unproject(pt: Point3) -> bool { - let proj = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0); +mod quickcheck_tests { + use na::{Point3, Perspective3, Orthographic3}; - let projected = proj.project_point(&pt); - let unprojected = proj.unproject_point(&projected); + quickcheck!{ + fn perspective_project_unproject(pt: Point3) -> 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) -> bool { - 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); + fn orthographic_project_unproject(pt: Point3) -> bool { + 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) + } } } diff --git a/tests/geometry/rotation.rs b/tests/geometry/rotation.rs index a6e6daf6..a2e91dec 100644 --- a/tests/geometry/rotation.rs +++ b/tests/geometry/rotation.rs @@ -1,6 +1,4 @@ -use std::f64; -use alga::general::Real; -use na::{self, Vector2, Vector3, Rotation2, Rotation3, Unit}; +use na::{Vector2, Vector3}; #[test] fn angle_2() { @@ -19,187 +17,193 @@ fn angle_3() { } #[cfg(feature = "arbitrary")] -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); +mod quickcheck_tests { + use std::f64; + use alga::general::Real; + use na::{self, Vector2, Vector3, Rotation2, Rotation3, Unit}; - 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. - pitch[(1, 1)] == 1.0 && // rotation wrt. y axis. - yaw[(2, 2)] == 1.0 && // rotation wrt. z axis. - yaw * pitch * roll == rpy - } + let rpy = Rotation3::from_euler_angles(r, p, y); - 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 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) -> 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) -> 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, b: Vector2) -> bool { - a.angle(&b) == b.angle(&a) - } - - fn angle_is_commutative_3(a: Vector3, b: Vector3) -> bool { - a.angle(&b) == b.angle(&a) - } - - /* - * - * Rotation matrix between vectors. - * - */ - fn rotation_between_is_anticommutative_2(a: Vector2, b: Vector2) -> 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, b: Vector3) -> 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) + roll[(0, 0)] == 1.0 && // rotation wrt. x axis. + pitch[(1, 1)] == 1.0 && // rotation wrt. y axis. + yaw[(2, 2)] == 1.0 && // rotation wrt. z axis. + yaw * pitch * roll == rpy } - 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, v3: Vector3) -> 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, b: Vector2) -> 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) + 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) } - else { - true + + /* + * + * Inversion is transposition. + * + */ + fn rotation_inv_3(a: Rotation3) -> 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, b: Vector3) -> 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) + fn rotation_inv_2(a: Rotation2) -> 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) } - else { - true + + /* + * + * Angle between vectors. + * + */ + fn angle_is_commutative_2(a: Vector2, b: Vector2) -> bool { + a.angle(&b) == b.angle(&a) + } + + fn angle_is_commutative_3(a: Vector3, b: Vector3) -> bool { + a.angle(&b) == b.angle(&a) + } + + /* + * + * Rotation matrix between vectors. + * + */ + fn rotation_between_is_anticommutative_2(a: Vector2, b: Vector2) -> 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, b: Vector3) -> 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, v3: Vector3) -> 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, b: Vector2) -> 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, b: Vector3) -> 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. - * - */ - fn new_rotation_2(angle: f64) -> bool { - let r = Rotation2::new(angle); + /* + * + * Rotation construction. + * + */ + fn new_rotation_2(angle: f64) -> bool { + 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) -> 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)) + relative_eq!(r.angle(), angle, epsilon = 1.0e-7) } - else { - r == Rotation3::identity() + + fn new_rotation_3(axisangle: Vector3) -> 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. - * - */ - fn powf_rotation_2(angle: f64, pow: f64) -> bool { - let r = Rotation2::new(angle).powf(pow); + /* + * + * Rotation pow. + * + */ + fn powf_rotation_2(angle: f64, pow: f64) -> bool { + let r = Rotation2::new(angle).powf(pow); - 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, 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 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)) + relative_eq!(r.angle(), pangle, epsilon = 1.0e-7) } - else { - r == Rotation3::identity() + + fn powf_rotation_3(axisangle: Vector3, 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() + } } } -); +} diff --git a/tests/geometry/similarity.rs b/tests/geometry/similarity.rs index 2364e992..feb8ba1b 100644 --- a/tests/geometry/similarity.rs +++ b/tests/geometry/similarity.rs @@ -1,9 +1,9 @@ +#![cfg(feature = "arbitrary")] #![allow(non_snake_case)] use alga::linear::{Transformation, ProjectiveTransformation}; use na::{Vector3, Point3, Similarity3, Translation3, Isometry3, UnitQuaternion}; -#[cfg(feature = "arbitrary")] quickcheck!( fn inverse_is_identity(i: Similarity3, p: Point3, v: Vector3) -> bool { let ii = i.inverse(); diff --git a/tests/geometry/unit_complex.rs b/tests/geometry/unit_complex.rs index 420b03c1..65837d8b 100644 --- a/tests/geometry/unit_complex.rs +++ b/tests/geometry/unit_complex.rs @@ -1,9 +1,8 @@ +#![cfg(feature = "arbitrary")] #![allow(non_snake_case)] use na::{Unit, UnitComplex, Vector2, Point2, Rotation2}; - -#[cfg(feature = "arbitrary")] quickcheck!( /* diff --git a/tests/linalg/balancing.rs b/tests/linalg/balancing.rs index 97af9ab6..73fd8230 100644 --- a/tests/linalg/balancing.rs +++ b/tests/linalg/balancing.rs @@ -1,9 +1,10 @@ +#![cfg(feature = "arbitrary")] + use std::cmp; use na::{DMatrix, Matrix4}; use na::balancing; -#[cfg(feature = "arbitrary")] quickcheck! { fn balancing_parlett_reinsch(n: usize) -> bool { let n = cmp::min(n, 10); diff --git a/tests/linalg/bidiagonal.rs b/tests/linalg/bidiagonal.rs index 30eddd4c..a96ce47c 100644 --- a/tests/linalg/bidiagonal.rs +++ b/tests/linalg/bidiagonal.rs @@ -1,7 +1,7 @@ +#![cfg(feature = "arbitrary")] + use na::{DMatrix, Matrix2, Matrix4, Matrix5x3, Matrix3x5}; - -#[cfg(feature = "arbitrary")] quickcheck! { fn bidiagonal(m: DMatrix) -> bool { if m.len() == 0 { diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs index 27c9b583..ceb44bdf 100644 --- a/tests/linalg/eigen.rs +++ b/tests/linalg/eigen.rs @@ -1,58 +1,61 @@ -use std::cmp; - -use na::{DMatrix, Matrix2, Matrix3, Matrix4}; +use na::DMatrix; #[cfg(feature = "arbitrary")] -quickcheck! { - fn symmetric_eigen(n: usize) -> bool { - let n = cmp::max(1, cmp::min(n, 10)); - let m = DMatrix::::new_random(n, n); - let eig = m.clone().symmetric_eigen(); - let recomp = eig.recompose(); +mod quickcheck_tests { + use std::cmp; + use na::{DMatrix, Matrix2, Matrix3, Matrix4}; - 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::::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 { - let n = cmp::max(1, cmp::min(n, 10)); - let mut m = DMatrix::::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_singular(n: usize) -> bool { + let n = cmp::max(1, cmp::min(n, 10)); + let mut m = DMatrix::::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) -> 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_4x4(m: Matrix4) -> 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) -> 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) -> 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) -> 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) -> 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) + } } } diff --git a/tests/linalg/full_piv_lu.rs b/tests/linalg/full_piv_lu.rs index 2cfa2c29..0f4b0a72 100644 --- a/tests/linalg/full_piv_lu.rs +++ b/tests/linalg/full_piv_lu.rs @@ -1,7 +1,4 @@ -use std::cmp; -use na::{DMatrix, Matrix3, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, - DVector, Vector4}; - +use na::Matrix3; #[test] fn full_piv_lu_simple() { @@ -42,114 +39,119 @@ fn full_piv_lu_simple_with_pivot() { } #[cfg(feature = "arbitrary")] -quickcheck! { - fn full_piv_lu(m: DMatrix) -> bool { - let mut m = m; - if m.len() == 0 { - m = DMatrix::new_random(1, 1); - } +mod quickcheck_tests { + use std::cmp; + use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, DVector, Vector4}; - let lu = m.clone().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_3_5(m: Matrix3x5) -> 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) -> 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) -> 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::::new_random(n, n); + quickcheck! { + fn full_piv_lu(m: DMatrix) -> bool { + let mut m = m; + if m.len() == 0 { + m = DMatrix::new_random(1, 1); + } let lu = m.clone().full_piv_lu(); - let b1 = DVector::new_random(n); - let b2 = DMatrix::new_random(n, nb); + 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 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)) + relative_eq!(m, lu, epsilon = 1.0e-7) } - return true; - } + fn full_piv_lu_static_3_5(m: Matrix3x5) -> 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) -> bool { - let lu = m.full_piv_lu(); - let b1 = Vector4::new_random(); - let b2 = Matrix4x3::new_random(); + relative_eq!(m, lu, epsilon = 1.0e-7) + } - let sol1 = lu.solve(&b1); - let sol2 = lu.solve(&b2); + fn full_piv_lu_static_5_3(m: Matrix5x3) -> 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); - 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)) - } + relative_eq!(m, lu, epsilon = 1.0e-7) + } - 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::::new_random(n, n); + fn full_piv_lu_static_square(m: Matrix4) -> 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); - let mut l = m.lower_triangle(); - let mut u = m.upper_triangle(); + relative_eq!(m, lu, epsilon = 1.0e-7) + } - // Ensure the matrix is well conditioned for inversion. - l.fill_diagonal(1.0); - u.fill_diagonal(1.0); - let m = l * u; + 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::::new_random(n, n); - let m1 = m.clone().full_piv_lu().try_inverse().unwrap(); - let id1 = &m * &m1; - let id2 = &m1 * &m; + let lu = m.clone().full_piv_lu(); + let b1 = DVector::new_random(n); + 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) -> bool { - let lu = m.full_piv_lu(); + 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)) + } - if let Some(m1) = lu.try_inverse() { + return true; + } + + fn full_piv_lu_solve_static(m: Matrix4) -> 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::::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 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) -> 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 + } } } } diff --git a/tests/linalg/hessenberg.rs b/tests/linalg/hessenberg.rs index 6fce8579..67206487 100644 --- a/tests/linalg/hessenberg.rs +++ b/tests/linalg/hessenberg.rs @@ -1,7 +1,8 @@ +#![cfg(feature = "arbitrary")] + use std::cmp; use na::{DMatrix, Matrix2, Matrix4}; - #[test] fn hessenberg_simple() { 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)) } - -#[cfg(feature = "arbitrary")] quickcheck! { fn hessenberg(n: usize) -> bool { let n = cmp::max(1, cmp::min(n, 50)); diff --git a/tests/linalg/lu.rs b/tests/linalg/lu.rs index 2e8cc0ce..7b823968 100644 --- a/tests/linalg/lu.rs +++ b/tests/linalg/lu.rs @@ -1,7 +1,4 @@ -use std::cmp; -use na::{DMatrix, Matrix3, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, - DVector, Vector4}; - +use na::Matrix3; #[test] fn lu_simple() { @@ -40,110 +37,115 @@ fn lu_simple_with_pivot() { } #[cfg(feature = "arbitrary")] -quickcheck! { - fn lu(m: DMatrix) -> bool { - let mut m = m; - if m.len() == 0 { - m = DMatrix::new_random(1, 1); - } +mod quickcheck_tests { + use std::cmp; + use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, DVector, Vector4}; - let lu = m.clone().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_3_5(m: Matrix3x5) -> 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) -> 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) -> 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::::new_random(n, n); + quickcheck! { + fn lu(m: DMatrix) -> bool { + let mut m = m; + if m.len() == 0 { + m = DMatrix::new_random(1, 1); + } let lu = m.clone().lu(); - let b1 = DVector::new_random(n); - let b2 = DMatrix::new_random(n, nb); + let (p, l, u) = lu.unpack(); + let mut lu = l * u; + p.inv_permute_rows(&mut lu); - 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)) + relative_eq!(m, lu, epsilon = 1.0e-7) } - return true; - } + fn lu_static_3_5(m: Matrix3x5) -> 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) -> bool { - let lu = m.lu(); - let b1 = Vector4::new_random(); - let b2 = Matrix4x3::new_random(); + relative_eq!(m, lu, epsilon = 1.0e-7) + } - let sol1 = lu.solve(&b1); - let sol2 = lu.solve(&b2); + fn lu_static_5_3(m: Matrix5x3) -> bool { + 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)) && - (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6)) - } + relative_eq!(m, lu, epsilon = 1.0e-7) + } - 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::::new_random(n, n); + fn lu_static_square(m: Matrix4) -> bool { + let lu = m.lu(); + let (p, l, u) = lu.unpack(); + let mut lu = l * u; + p.inv_permute_rows(&mut lu); - let mut l = m.lower_triangle(); - let mut u = m.upper_triangle(); + relative_eq!(m, lu, epsilon = 1.0e-7) + } - // Ensure the matrix is well conditioned for inversion. - l.fill_diagonal(1.0); - u.fill_diagonal(1.0); - let m = l * u; + 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::::new_random(n, n); - let m1 = m.clone().lu().try_inverse().unwrap(); - let id1 = &m * &m1; - let id2 = &m1 * &m; + let lu = m.clone().lu(); + let b1 = DVector::new_random(n); + 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) -> bool { - let lu = m.lu(); + 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)) + } - if let Some(m1) = lu.try_inverse() { + return true; + } + + fn lu_solve_static(m: Matrix4) -> 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::::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 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) -> 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 + } } } } diff --git a/tests/linalg/qr.rs b/tests/linalg/qr.rs index 8bdc3ce8..394a8b42 100644 --- a/tests/linalg/qr.rs +++ b/tests/linalg/qr.rs @@ -1,8 +1,9 @@ +#![cfg(feature = "arbitrary")] + use std::cmp; use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, DVector, Vector4}; -#[cfg(feature = "arbitrary")] quickcheck! { fn qr(m: DMatrix) -> bool { let qr = m.clone().qr(); diff --git a/tests/linalg/real_schur.rs b/tests/linalg/real_schur.rs index 3ef32297..40f8a011 100644 --- a/tests/linalg/real_schur.rs +++ b/tests/linalg/real_schur.rs @@ -1,5 +1,4 @@ -use std::cmp; -use na::{DMatrix, Matrix2, Matrix3, Matrix4}; +use na::{DMatrix, Matrix3, Matrix4}; #[test] @@ -15,49 +14,54 @@ fn schur_simpl_mat3() { } #[cfg(feature = "arbitrary")] -quickcheck! { - fn schur(n: usize) -> bool { - let n = cmp::max(1, cmp::min(n, 10)); - let m = DMatrix::::new_random(n, n); +mod quickcheck_tests { + use std::cmp; + use na::{DMatrix, Matrix2, Matrix3, Matrix4}; - 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::::new_random(n, n); - if !relative_eq!(&vecs * &vals * vecs.transpose(), m, epsilon = 1.0e-7) { - println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose()); + let (vecs, vals) = m.clone().real_schur().unpack(); + + 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) -> bool { + let (vecs, vals) = m.clone().real_schur().unpack(); - fn schur_static_mat2(m: Matrix2) -> bool { - let (vecs, vals) = m.clone().real_schur().unpack(); - - let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7); - if !ok { - println!("{:.5}{:.5}", vecs, vals); - println!("Reconstruction:{}{}", m, &vecs * &vals * vecs.transpose()); + let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7); + if !ok { + println!("{:.5}{:.5}", vecs, vals); + println!("Reconstruction:{}{}", m, &vecs * &vals * vecs.transpose()); + } + ok } - ok - } - fn schur_static_mat3(m: Matrix3) -> bool { - let (vecs, vals) = m.clone().real_schur().unpack(); + fn schur_static_mat3(m: Matrix3) -> bool { + let (vecs, vals) = m.clone().real_schur().unpack(); - let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7); - if !ok { - println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose()); + let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7); + if !ok { + println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose()); + } + ok } - ok - } - fn schur_static_mat4(m: Matrix4) -> bool { - let (vecs, vals) = m.clone().real_schur().unpack(); + fn schur_static_mat4(m: Matrix4) -> bool { + let (vecs, vals) = m.clone().real_schur().unpack(); - let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7); - if !ok { - println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose()); + let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7); + if !ok { + println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose()); + } + ok } - ok } } diff --git a/tests/linalg/solve.rs b/tests/linalg/solve.rs index 3940c287..f960fb0a 100644 --- a/tests/linalg/solve.rs +++ b/tests/linalg/solve.rs @@ -1,3 +1,5 @@ +#![cfg(feature = "arbitrary")] + use na::{Matrix4, Matrix4x5}; fn unzero_diagonal(a: &mut Matrix4) { @@ -8,7 +10,6 @@ fn unzero_diagonal(a: &mut Matrix4) { } } -#[cfg(feature = "arbitrary")] quickcheck! { fn solve_lower_triangular(a: Matrix4, b: Matrix4x5) -> bool { let mut a = a; diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index df107694..5b09c7c8 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -1,142 +1,143 @@ -use std::cmp; - -use na::{DMatrix, Matrix2, Matrix3, Matrix4, Matrix6, Matrix5x2, Matrix5x3, Matrix2x5, Matrix3x5, - DVector}; - +use na::{DMatrix, Matrix6}; #[cfg(feature = "arbitrary")] -quickcheck! { - fn svd(m: DMatrix) -> 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); +mod quickcheck_tests { + use std::cmp; + use na::{DMatrix, Matrix2, Matrix3, Matrix4, Matrix5x2, Matrix5x3, Matrix2x5, Matrix3x5, DVector}; - println!("{}{}", &m, &u * &ds * &v_t); + quickcheck! { + fn svd(m: DMatrix) -> 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) && - relative_eq!(&u * ds * &v_t, recomp_m, epsilon = 1.0e-5) && - relative_eq!(m, recomp_m, epsilon = 1.0e-5) - } - else { - true - } - } + println!("{}{}", &m, &u * &ds * &v_t); - fn svd_static_5_3(m: Matrix5x3) -> 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) -> 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) -> 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) -> 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) -> 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) -> 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) -> 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) + s.iter().all(|e| *e >= 0.0) && + relative_eq!(&u * ds * &v_t, recomp_m, epsilon = 1.0e-5) && + relative_eq!(m, recomp_m, epsilon = 1.0e-5) } else { - println!("{}", &m * &pinv); - (m * pinv).is_identity(1.0e-5) + true } } - else { + + fn svd_static_5_3(m: Matrix5x3) -> 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) -> 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) -> 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) -> 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) -> 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) -> 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) -> 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::::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 } } - - 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::::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. diff --git a/tests/linalg/tridiagonal.rs b/tests/linalg/tridiagonal.rs index 2128d983..6db86aef 100644 --- a/tests/linalg/tridiagonal.rs +++ b/tests/linalg/tridiagonal.rs @@ -1,9 +1,9 @@ +#![cfg(feature = "arbitrary")] + use std::cmp; use na::{DMatrix, Matrix2, Matrix4}; - -#[cfg(feature = "arbitrary")] quickcheck! { fn symm_tridiagonal(n: usize) -> bool { let n = cmp::max(1, cmp::min(n, 50));