From a52b07957845884da10704b94ab68e0358c32a69 Mon Sep 17 00:00:00 2001 From: Andreas Longva Date: Thu, 27 Apr 2017 22:19:32 +0200 Subject: [PATCH] Relax invertibility test in try_inverse() The previous implementation of try_inverse() used an approximate check of the determinant against 0 for small matrices to determine if the matrix was invertible. This is not a reliable test, and may fail for perfectly invertible matrices. This change simply makes the test criterion an exact comparison instead. --- src/core/inverse.rs | 9 ++++--- tests/matrix_inverse.rs | 58 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 4 deletions(-) diff --git a/src/core/inverse.rs b/src/core/inverse.rs index ed066e49..488c5b5b 100644 --- a/src/core/inverse.rs +++ b/src/core/inverse.rs @@ -45,18 +45,19 @@ impl SquareMatrix match dim { 0 => true, 1 => { - if relative_eq!(self.get_unchecked(0, 0), &N::zero()) { + let determinant = self.get_unchecked(0, 0).clone(); + if determinant == N::zero() { false } else { - *self.get_unchecked_mut(0, 0) = N::one() / self.determinant(); + *self.get_unchecked_mut(0, 0) = N::one() / determinant; true } }, 2 => { let determinant = self.determinant(); - if relative_eq!(&determinant, &N::zero()) { + if determinant == N::zero() { false } else { @@ -94,7 +95,7 @@ impl SquareMatrix m12 * minor_m11_m23 + m13 * minor_m11_m22; - if relative_eq!(&determinant, &N::zero()) { + if determinant == N::zero() { false } else { diff --git a/tests/matrix_inverse.rs b/tests/matrix_inverse.rs index a37b54ff..fe1af112 100644 --- a/tests/matrix_inverse.rs +++ b/tests/matrix_inverse.rs @@ -61,3 +61,61 @@ fn matrix5_try_inverse() { assert_relative_eq!(a_inv, expected_inverse, max_relative=1e-4); } + +#[test] +fn matrix1_try_inverse_scaled_identity() { + // A perfectly invertible matrix with + // very small coefficients + let a = Matrix1::new(1.0e-20); + let expected_inverse = Matrix1::new(1.0e20); + let a_inv = a.try_inverse().expect("Matrix is invertible"); + + assert_relative_eq!(a_inv, expected_inverse); +} + +#[test] +fn matrix2_try_inverse_scaled_identity() { + // A perfectly invertible matrix with + // very small coefficients + let a = Matrix2::new(1.0e-20, 0.0, + 0.0, 1.0e-20); + let expected_inverse = Matrix2::new(1.0e20, 0.0, + 0.0, 1.0e20); + let a_inv = a.try_inverse().expect("Matrix is invertible"); + + assert_relative_eq!(a_inv, expected_inverse); +} + +#[test] +fn matrix3_try_inverse_scaled_identity() { + // A perfectly invertible matrix with + // very small coefficients + let a = Matrix3::new(1.0e-20, 0.0, 0.0, + 0.0, 1.0e-20, 0.0, + 0.0, 0.0, 1.0e-20); + let expected_inverse = Matrix3::new(1.0e20, 0.0, 0.0, + 0.0, 1.0e20, 0.0, + 0.0, 0.0, 1.0e20); + let a_inv = a.try_inverse().expect("Matrix is invertible"); + + assert_relative_eq!(a_inv, expected_inverse); +} + +#[test] +fn matrix5_try_inverse_scaled_identity() { + // A perfectly invertible matrix with + // very small coefficients + let a = Matrix5::new(1.0e-20, 0.0, 0.0, 0.0, 0.0, + 0.0, 1.0e-20, 0.0, 0.0, 0.0, + 0.0, 0.0, 1.0e-20, 0.0, 0.0, + 0.0, 0.0, 0.0, 1.0e-20, 0.0, + 0.0, 0.0, 0.0, 0.0, 1.0e-20); + let expected_inverse = Matrix5::new(1.0e+20, 0.0, 0.0, 0.0, 0.0, + 0.0, 1.0e+20, 0.0, 0.0, 0.0, + 0.0, 0.0, 1.0e+20, 0.0, 0.0, + 0.0, 0.0, 0.0, 1.0e+20, 0.0, + 0.0, 0.0, 0.0, 0.0, 1.0e+20);; + let a_inv = a.try_inverse().expect("Matrix is invertible"); + + assert_relative_eq!(a_inv, expected_inverse); +}